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Abstract 



The thesis deals with the first stage of planet formation, namely dust coagulation from 
micron to millimeter sizes in circumstellar disks. For the first time, we collect and com- 
pile the recent laboratory experiments on dust aggregates into a collision model that 
can be implemented into dust coagulation models. Wc put this model into a Monte 
Carlo code that uses representative particles to simulate dust evolution. Simulations are 
performed using three different disk models in a local box (OD) located at 1 AU distance 
from the central star. We find that the dust evolution does not follow the previously 
assumed growth-fragmentation cycle, but growth is halted by bouncing before the frag- 
mentation regime is reached. We call this the bouncing barrier which is an additional 
obstacle during the already complex formation process of planetesimals. The absence of 
the growth-fragmentation cycle and the halted growth has two important consequences 
for planet formation. 1) It is observed that disk atmospheres are dusty throughout their 
lifetime. Previous models concluded that the small, continuously produced fragments 
can keep the disk atmospheres dusty. We however show that small fragments are not 
produced because bouncing prevents fragmentation. 2) As particles do not reach the 
fragmentation barrier, their sizes are smaller compared to the sizes reached in previous 
dust models. Forming planetesimals from such tiny aggregates is a challenging task. 
Wc decided to investigate point 1) in more detail. A vertical column of a disk (ID) is 
modeled including the sedimentation of the particles. We find that already intermediate 
levels of turbulence can prevent particles settling to the midplane. We also find that, 
due to bouncing, the particle size distribution is narrow and homogenous as a function 
of height in the disk. This finding has important implications for observations. If it is 
reasonable to assume that the turbulence is constant as a function of height, the particles 
observed at the disk atmospheres have the same properties as the ones at the midplane. 
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Zusammenfassung 



Diese Arbeit befasst sich mit der friihesten Phase der Planetenentstehung, namlich der 
Koagulation von mikrometer- hin zu millimetergrofien Staubpartikeln in zirkumstellaren 
Scheiben. Als erste Studie dieser Art simulieren wir die Staubentwicklung in 'representa- 
tive particle' Monte-CarloSimulationen unter Verwendung eines Kollisionsmodells, das 
die neuesten Laborexperimente beriicksichtigt. Die Simulationen verwenden drei ver- 
schicdcne Scheibenmodcllc in cincr lokalen Box (OD) in cincm Abstand von 1 AU vom 
Zentralstern. Unsere Ergebnisse zeigen, dass die Staubentwicklung nicht dcm bislang 
angenommenen Wachstums-Fragmentations-Zyklus folgt, sondcrn dass das Wachstum 
von abprallenden Stofien aufgehalten wird, bevor es das Fragnicntationsrcgime erreicht. 
Wir bezeichnen dies als 'bouncing barrier', cin weiteres Hindcrnis im ohnehin schon 
komplexen Entstehungsprozess von Planetesimalen. Die Abwesenheit des Wachstums- 
Pragmentations-Zyklus und das unterbundene Teilchenwachstum haben zwei wichtige 
Konsequenzen fiir die Entstehung von Planeten: 1) Beobachtungen zeigen, dass die 
Atmospharen von Scheiben wahrend ihrer gesamten Lebenszeit staubig sind. Bish- 
erige Modelle folgerten dass kontinuierliche Fragmentation diese kleinen Staubteilchen 
produziert und dadurch die Scheibenatmosphare "staubig" halt. Unsere Ergebnisse 
zeigen jedoch, dass kleine Fragmente gar nicht erst produziert werden, well die Frag- 
mentationsgrenze nicht erreicht wird. 2) Da Teilchen die Fragmcntationsbarriere nicht 
erreichcn, blcibcn sic klcincr als in bisherigen Modellen. Die Entstehung von Plan- 
etesimalen aus solch kleinen Staubaggregaten ist eine hcrausfordcrungsvollc Aufgabc. 
Wir haben uns mit Punkt 1) naher befasst. Hierzu modcUicrcn wir cincn vcrtikalen 
Schnitt (ID) durch die Scheibe unter Berticksichtigung von Staubsedimentation. Unsere 
Ergebnisse zeigen, dass schon eine moderat ausgeprgte Turbulenz die Sedimentation zur 
Mittelebene unterbinden kann. Des Weiteren fanden wir her aus, dass die Verteilung 
der TeilchengroBe schmal und eine homogene Funktion der Hohe iiber der Mittelebene 
ist. Dies hat wichtige Auswirkungen fiir Beobachtungen: Unter der Annahme, dass 
die Turbulenz hohenunabhangig ist, haben die in der Scheibenatmosphare beobachteten 
Teilchen dieselben Eigenschaften wie diejenigen in der Mittelebene. 
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Chapter 1 



Introduction 



Perhaps the most important astronomical discovery of the second half of the last century 
was the detection of planets orbiting around nearby sun-like stars. This discovery and 
the since then vastly growing number of exoplanet discoveries further trigger our interest 
in their formation, their atmospheres, chemical and geophysical evolution, etc. 

The theory of planet formation needs to successfully explain all the different types of 
exoplanets ranging from the Hot-Jupiters, which revolve around their star in a period of 
days, the presumably rocky Super-Earth planets, and giant planets orbiting more than 
100 AU distance from the central star. This variety and the complex processes of planet 
formation make the field fascinating. 

The formation of stars starts with the gravitational collapse of a dense molecular cloud. 
Due to the initial angular momentum of the cloud, most of the infalling matter will not 
fall directly onto the protostar, but it forms a disk around it. The matter in the dense 
cloud is twofold: roughly 99% (in mass) of the matter is present in the form of gases, 
the rest consists of solid, sub-micron sized dust particles. This solid matter serves as 
the building blocks of planets. We see that the way from the sub-micron sized particles 
(10~^^ g) to a planet of several lOOO's of kilometers in radius (10^^ g) is a long one (see 
Sec. 11.11 for more details). We do not attempt to go the whole way in this thesis, but 
rather concentrate on the initial stages of growth until millimeter-centimeter (1 g) in 
size. 

To understand planet formation and dust evolution, first we need to know the typical 
physical conditions of the environment where it happens, i.e., the properties of proto- 
planetary disks. Our solar system provides only vague clues regarding its formation ^ 
history, as the solar nebula has long dissipated: the distribution and properties of the 
planets, minor bodies, satellites and meteorites all contribute to our understanding. 
These results are reviewed in Sec. 11.2.11 Observations of existing protoplanetary disks 

*The images at the lower right corners of this thesis show the snapshots of the sedimentation simu- 
lation presented in Chapter 15.3.51 Fig. 15.81 as a flip-cartoon. 
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in all wavelengths from UV to millimeter revealed the basic properties of planet-forming 
disks such as their typical masses, lifetimes, physical conditions (density and tempera- 
ture structure), important processes such as accretion, photoevaporation, etc. We review 
these processes in Sec. 11.2.21 

Once the typical physical properties of the disk are known, we need to understand the 
driving force of particle growth. The solid particles in the gas disk feel the drag from the 
gas, and have relative velocities. As a rule of thumb, we can say that the relative velocity 
is increasing with mass until the particles reach a meter in size. Due to the (random 
and systematic) relative velocities, the particles can collide and grow, therefore it is 
crucial to understand what sources of relative velocity are there. One should also keep 
in mind that low velocity collisions produce self-similar so called fractal structures, while 
intermediate and high velocities lead to deformation (restructure and/or destroy) of the 
aggregates. Such deformation changes the aerodynamical properties of the particles and 
thus in turn the relative velocity (see Chapters 14.2.21 and I5.2.3P . 

We have to know the outcome of each collision: will the aggregates stick, bounce or 
fragment? During the last fifteen years, a huge amount of laboratory data was collected 
about the collisions of dust aggregates. We know that the monomers (sub-micron sized 
particles) stick together due to surface forces (van der Waals attraction for silicates and 
molecular dipole interaction for ices). We can measure the strength of the surface forces, 
furthermore it is possible to produce fractal aggregates and perform experiments with 
porous, but non- fractal structures in a wide velocity range (from mm/s to several lO's 
of m/s, see Sec. II. 3p . 

The laboratory experiments are crucial for our understanding of the microphysics of 
collisions as we cannot gain information about it in any other way. The 'molecular 
dynamics' models can also be used to simulate collisions of aggregates, however these 
models also rely on measured quantities, such as surface energy. Young's modulus, etc. 
The difficulty with the experiments is that they produce results which are not always 
easily usable in theoretical modeling, as the results are rather complex. One of the main 
goals of this thesis is to collect all the available laboratory data and compile it into a 
collision model (Chapter [3]). 

Once we know what happens during individual collisions, we have to calculate how the 
entire population of dust particles evolves. Traditional methods integrate the Smolu- 
chowski equation with only one particle property: the particle mass (the equation, and 
such solvers are introduced and described in Chapter 12. 4p . Although such methods are 
fast, thus the whole disk can be modeled, they have difficulty including any additional 
dust property. Therefore we use a Monte Carlo method (described in Chapter [2]), which 
is flexible and can follow the porosity of the aggregates as well (which is a dust-parameter 
in the experiment-based collision model), and it is straightforward to include a collision 
model with arbitrary complexity. The price we pay for this flexibility is the speed of 
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the code. We can only follow dust evolution in a local box (see Chapter d]) or in a ID 
vertical column (Chapter [5|) so far. 

One should also realize how entangled dust coagulation is with the other ongoing pro- 
cesses in disks. A truly self-consistent model should take into account the gas and dust 
evolution, how the surfaces of dust grains influence the chemistry and ionization fraction 
of the gas, how this ionization fraction regulates the coupling of the gas with the mag- 
netic fields, how this coupling determines the strength of turbulence, and finally closing 
the loop, how the strength of turbulence determines the dust evolution. The opacity of 
dust influences the observed properties of disks, the dust also sets the temperature of the 
gas at the main body of the disk (except at the low density upper layers), and the dust 
particles (if sufficiently concentrated) can influence the gas dynamics as well. Although 
such a complex model, which takes all these processes self-consistently into account, does 
not exist yet, the effects of the individual processes are at least approximately known 
(see Sec. fLl) . 



1.1 Planet formation in a nutshell 



Star formation mostly happens in star clusters in which roughly half of the stars are 
part of a binary or small multiple systems. The formation of an isolated star is easier 
to understand (no initial fragmentation of the parent cloud onto multiple objects, no 
environmental effects, such as stellar encounters and strong external radiation, have to 
be treated), therefore isolated star formation is better-studied and understood than star 
formation in an active environment. 

Star formation starts with the collapse of the parent cloud. This cloud rotates and 
has too much angular momentum to collapse directly onto the protostar, therefore a 
disk forms around the cen tral object. It is obser ved that matter is accreted from the 



disk onto the central star (Hartmann et al 



19981 ). To acc rete, angular momentum ha s 



to be redistributed wit h in th e disk (viscous spreading - IShakura &: SunvaevI (|l973l ) ; 
Lynden-Bell Sz Pringlel (1974)), o r angular mornenturn has to be lost (e.g., by disk 
winds - Blandford &: Payn^ ( 1982 ): Lubow et al. ( 1994 )). These processes happen on 
a much longer timescale than the orbital period, therefore in many applications it is 
a reasonable assumption that disks are in quasi-equilibrium. How angular momentum 
is lost /redistributed is one of the most actively studied area of star formation and the 
physics of disks (see Sec. 11.41 for more details). 



The parent cloud also contains dust particles of sub-micrometer in size (jMathis et al. 



19771 ). It is usually assumed that the dust-to-gas ratio in the interstellar matter is 1:100 
by mass. This ratio is also inherited by the disk and this solid material (iron, silicates 
and ices) is the ingredient of terrestrial planets and the cores of some giant planets. 
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Curre ntly there are two leading giant planet formation models: disk instability (jBossl . 
19971 ) and core accretion. These models are often thought of being competitive to each 
other, but one can imagine that these mechanisms operate at different parts of the disk 
simultaneously. 



The core accretion model (iMizuno . 



198C 



Pollack et al. 



19961 ) starts with the assembly 



of the protoplanetary core. This process converts the sub-micron sized dust particles 
into a core which has some thousands of kilometers in size. The assembly of the core 
happens in three steps. 1.) The dust particles grow by two-body interactions where the 
gravitational interaction between the bodies is negligible, in other words dust coagula- 
tion. Coagulation can easily produce aggregates of millimeter in size. 2.) How these 
aggregates are exactly converted into planetesimals is one of the key problems in planet 
formation. This step bears from many problems (radial drift barrier, fragmentation bar- 
rier, bouncing barrier), which are discussed in detail in the upcoming chapters of this 
thesis. Planetesimals could form via co agulation if the envir onment is quiet, meaning 



that turbulence and radial drift are low (Brauer et al 



2008bl ). Or millimeter-decimeter 



sized particles could be concentrated in turbulent eddies and/or in vortices, where these 
concentrations are further enhanced by self-gravity of the particles. These particles could 
directly coalesce into kilometer sized planetesimals via many-body interactions (see the 
discussion in Chapter 14.5. 3p . Such a process cannot be modeled by the methods of this 
thesis, as we can only follow two-body interactions. 3.) Once planetesimals form, these 
planetesimals grow further due to gravitational agglomeration. 

Once the core is formed, and gas is still present in the disk, an atmosphere is gathered 
around the core. Initially the envelope around the core is in hydrostatic equilibrium. 
The energy gained by the impacting planetesimals and the gravitational potential energy 
released due to the contraction of the atmosphere is in equilibrium with the energy loss 
due to convection and radiative diffusion. When, however, a critical mass is exceeded, 
runaway gas accretion starts. This runaway accretion ends when no gas is present around 
the planet due to the presence of a gap or disk dissipation. After this stage, the planet 
goes through a Kelvin-Helmholtz contraction. 

The timescales involved are a se rious uncertainty of t he co re accretion model. A typic al 



disk lifetime is ~ 10^ yrs (e.g., iHaisch et al.l (|200ll ) and ISicilia-Aguilar et al.l (|2006l )) 
during which the core has to be assembled and sufficient amount of gas has to be accreted. 



200 



tI) and interact gravitationally due to resonances ( Lecar et al. . 2001 ). Once 



Once the planets are formed, thev can migrate in the gas disk (iMasseti. 120081: Pap aloizou 
et al., 



the gas disk is dispersed e.g.. by photo-evaporation fdue to external radiation - Adams 



et al. ( 2004). and due to radiation from the central star - [Alexander et al.l (120061): Gorti 
et al. (|2009l )). a so-called debris disk can still be present around the young star. 
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1.2 Observations of planet-forming systems 

The planetary system that we can examine in the greatest detail is our solar system. We 
know the positions and physical properties of the eight planets, we are able to observe 
the asteroid belt, minor bodies, and satellites in the solar system, furthermore we can 
directly examine meteorites in the laboratory, determine their absolute and relative ages. 
All this information contains hints about the formation history of the solar system. 

We can also observe other planet-forming systems in various wavelengths ranging from 
millimeter to UV, and we can use molecular spectral lines to trace the different molecular 
species. These observations put constraints on the physical properties of disks. 



1.2.1 The solar system 



We do not know the exact mass of the dust and gas from which the planets in the solar 
system formed. However, it is possible to derive a lower mass limit using the masses and 
positions of the solar system planets. This nebula model is called the minimum mass 
solar nebula model (MMSN) introduced by IWeidenschillingI (|l977l ). In the model, the 
masses of heavy elements of the planets are mixed with hydrogen and helium to reach 
solar compositions. The solar system is divided into concentric annuli each centered 
around the planet and extending half way until the next planet. The matter is then 
spread homogeneously in these annuli to obtain the gas surface density at the location 
of each planet. Performing these steps we get that the surface de nsity scales as S(r) tx 
^-3/2_ rjij^g most commonly used normalization is ( Havashil . 1 1981 ) 



S(r) = 1.7 X 10^ 



1 AU 



-3/2 



g cm 



(1.1) 



This is a lower limit of the solar nebula because it assumes that all the solids presented 
in the disk were incorporated into the planets. The model also assumes that the planets 
were formed in their current location, which might not be true due to the effects of mi- 
gration. There are attempts to u pdate this mo del by taking into account the migration 
of planets. Such calculations bv iDeschI (j2007l ) found that S(r) oc r~^-^, although this 
model is debate d in the literature. Hqw'ever, based on millimeter observations of proto- 
planetary disks, [Andrews &: Williamsl ( 2007 ) found that S(r) oc r~^/^. Their findings 
are probably more representative at the outer parts of the disks (~ 100 AU) and not at 
the inner parts of the disks. We conclude that the mass distribution in the early solar ^v.^'^^p ssom -lujou 
system is still very much debated. 



p IT) 

d d 



The regular planetary satellites of Jupiter (the Galilean satellites) have similar masses, 
tight prograde orbits which lie close to the equatorial plane of the planet, and they 
are trapped in mean motion resonances. This suggests that these satellites formed in a 
disk which surrounded the planet similarly to the primordial solar nebula surrounding 
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the early sun. On the other hand, Saturn has only one big regular satellite, Titan. The 
differences of the Jovian and Saturnian systems can be explained by the quick truncation 
of infalling matter in the Jovian system, which is caused by Jupiter opening a gap in the 
solar nebula. Due to the lower mass of Saturn, the Saturnian disk had a longer lifetime, 
therefore one single big moon could assemble (j Sasaki et al.l . 12010 ) . 



Other irregular satellites orbit in a larger distance from the host planet, sometimes in 
retrograde orbits. Such satellites were probably captured by the planet. 

The solar system also contains many minor bodies. In the inner solar system the asteroid 
belt is prominent. The distribution of the semi major axes of the asteroid belt objects is 
not homogenous, instead, gaps can be observed (Kirkwood gaps) and the asteroid belt 
can be divided into three regions: inner belt (distance < 2.5 AU - 3:1 resonance with 
Jupiter), central belt (distance between 2.5 and 3.81 AU - 2:1 resonance), and outer 



belt (beyond 3.81 AU). These asteroid s were heav i. 



they could not assemble into a planet (iPetit et al 



y pe rturbed by Jupiter therefore 



20011). Their m ass is also greatly 



depleted compared to the primordial mass (jWeidenschillingl . 119771 ). and the different 
asteroid types are radially mixed. It seems that in order to explain all these properties 
of the asteroid belt, a combination of sweeping secular resonances from the migrating 
Jupiter and Saturn, and ernbedded planetary embryos are needed that excite and scatter 
one another ([O'Brien et al.l . 120071 ) . 



Perhaps more interesting in the context of the primordial solar nebula are the Kuip er 
Belt Objects in the outer solar system (KBO - for a review see lLuu &: JewittI (120021 )). 
These objects have several classes like the classical KBOs (with low eccentricity), res- 
onant KBOs (in 4:3, 3:2, 2:1 resonance with Neptune ~ intermediate eccentricity), and 
the scattered KBOs (with eccentricities around 0.5). As we see, the Kuiper belt is dy- 
namically excited by Neptune. Two other properties of the Kuiper belt are important. 
1.) The Kuiper belt has a mass of 0.1 Earth mass, which is surprisingly low. Accretion 
models predict that a mass of 1 Earth mass must hav e existed to explain the growth of 
the obj ects we see now (see e.g. iKenyon &: Luul (jl998l )). 2.) The Kuiper belts ends near 
50 AU. 



Gomes et al. 



(120041 ) examined two scenarios which could explain the low total 
mass of KBOs: the vast majority of KBOs crossed orbits with Neptune, therefore their 
orbits were scattered; fragmentation into dust removed most of the mass of the Kuiper 
belt. They concluded that none of these two scenarios are likely, instead they proposed 
that the protoplanetary disk possessed an edge about 30 AU. This edge is responsible for 
stopping the outward migration o f Neptune, and durin g this migration, Neptune could 
have pushed the KBOs outwards (jLevison et al. , 2008 ). 



We have the unique possibility to measure the ages of solar system rocks accurately 
with the means of radioactive dating. Although the accuracy of dating in the solar 
system is unmatched with any astronomical observations, one must keep in mind that 
this accuracy is achieved for a single (and rather special) system. One can calculate 
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absolute and relative ages of meteorites u s ing lo n g-lived and short-live d radionuclides 
respectively (for a review see 



Russell et al 



mm; 



Wadhwa et al, 



(120071)) 



Primitive meteorites (chondrites) were never differentiated, they are relatively unaltered 
since their time of formation, therefore they preserve the early history of the solar 
system. Most primitive meteorites contain small inclusions which were heated to high 
temperatures, during which the amorphous dust became crystallized. Spectral signatures 
of crystalline dust in protoplanetary disks are common, such crystalliza tion process was 



obser ved 'in action' after an outburst of a Sun-like young star, EXLupi (jAbraham et al, 
20od ). 



Chondrules are the most abundant type of inclusions, more than 70% of the volume 
of primitive meteorites are chondrules. CAIs (calcium aluminum inclusions) are rarer, 
but they were subject to a more extreme heating event. The formation, specifically the 
heati ng mechanism for chondrules and CAIs is the subject of active research f Connollv 



et al., I2OO6I ). The three main hypotheses are: 1.) heating near the young Sun (X-wind 
model) with strong outward transport; 2.) shock waves in the gas disk; 3.) collisions 
between planetesimals and/or protoplanets. 

CAIs are the oldest object s in the solar s y stem with 4567.11 it 0.16 Myr as determined 
from ^''^Pb-^'^^Pb dating ( Russell et al.l . I2OO6I ). This age pinpoints a specific event, 
namely the solidification of CAIs. One has to keep in mind that other events, like the 
collapse of the Sun's parent cloud happened earlier, but such events cannot be dated 
precisely. The absolute age determination only recently became accurate enough to 
reliably determine the time interval between the formation of CAIs and chondrules. 
Such measurements suggest that some Myr passed between the formation of the oldest 
CAIs and the formation of chondrules. Absolute age determination methods assume that 
the abundance of parent and daughter isotopes is only altered via radioactive decay. 

The relative ages between chondrules and CAIs can be determined using another method 
using short-lived radionuclides such as ^®A1. The method also assumes that the abun- 
dance of parent and daughter isotopes is altered only via radioactive decay. Furthermore, 
it is also assumed that ^^Al was uniformly distributed in the disk (isotopic equilibration), 
otherwise the differences in the orig inal 26A1/27A1 ratio can indicate different local for- 
mation environments. The ^^Al method also yields to a relative age of 1-2 Myr between 
the CAIs and chondrules. 

This age spread fits within the observed lifetime of the disk, and (as the chondrules and 
CAI's coexist in meteorites) it suggests that the formation of planetesimals were either 
delayed or ongoing for several Myr. A rapid plane t esima l formation in less than 1 Myr 
is not supported by meteoritic data (jRussell et al.l . 120061 ). 
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1.2.2 Observations of protoplanetary systems 



Most of the information about protoplanetary disks is obtained through infrared (IR) 
measurements. Disks have an IR excess due to the presence of dust particles. These 
particles can scatter or absorb the radiation of the central star, they also reemit thermal 
radiation in the IR. Near-IR (NIR) wavelengths map the warm dust close to the star 
(order of 1 AU or smaller), which originates from the upper layers of the disk. Far- 
IR (FIR) radiation maps colder dust further away from the star. The spectral energy 
distribution (SED) can be used to fit the disk parameters (such as disk geometry, dust 
composition, temperature an d density structure - for a 2D radiative transfer model see 
Dulleinond &: DominikI (l2004l'l. a review about disk modeling can be found in Dullemond 
et al. (|2007h ). 



The IR excess of disks can be used to determin e the typical lifetime of disks (e.g., 
(j200ll ) and lSicilia-Aguilar et al.l (j2006l )). One can measure the disk fraction 



Haisch et al 



of young stars (e.g., how many percentage of stars have IR excess) in different young 
clusters and correlate this disk fraction with the estimated age of the cluster. Following 
this procedure it was concluded that only half of the young stars have disks in a cluster 
of 2 Myrs age. 

A prominent feature of the IR spectra is the 10 micron feature, which originates from 
small silicate particles (order of sub-micron in size) at the upper layers of disks. The 
shape of the 10 micron feature c ontai ns valuable information about the composition of 
these grains (Ivan Boekel et al.l. EooJ, although some caution is required when inter- 



preting the data (jJuhasz et al. . l200i^) . The most interesting aspect of the 10 micron 
feature for t he topic of this thesis is that it can provide evidence for grain growth ( van 
Boekel et al.. 
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I" 



Some sources have a flatter and broader spectral feature which 
suggests that in these sources the particles grew to a few micron in size at the upper 
layers. Theoretical models including particle settling, coagulation, and radiative trans- 
fer cannot vet convincinglv reproduce this aspect of the 10 micron feature ( Dullemond 
& Dominik, 
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Sedimentation driven coagulation is the topic of Chapter where 
these physical processes are discussed in detail. 

Optical and NIR scattered light images contain lot of information about disks. The 
used observational technique differs for different disk inclinations: coronagraphic mea- 
surements are used for high and intermediate inclinations; observing optically thick 
edge-on disks require high spatial resolution, but not high contrast, therefore adaptive 
optics systems are advantageous; the Orion nebula provides a unique opportunity to 
observe silhouette disks (the disks appear dark due to the bright background of the neb- 
ula). The information that these images provide also depends on the inclination. From 
face-on disks we can determine the ellipticity of the disk, the dependence of the surface 
brightness on the radius, and the ratio of scattered and unscattered light (the relative 
brightness of the disk and the star). From intermediate inclinations (when the central 
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star is still visible), one can determined the relative brightness of the disk and the star, 
and the outer radius of the disk (if both the upper and lower parts of the nebulae are 
visible). In case of the edge on disks the inclination can be very precisely determined, 
and the effective scale-height of dust in the outer d isk. A detailed discussion of scattered 
light images can be found in 



Watson et al 



(120071 ). 



Millimeter observations are a useful tool to estimate the mass of the solid material in 
the disk and these observations provide evidence for grain growth at the outer disk. 
At millimeter and sub-millimeter wavelengths (assuming the Rayleigh-Jeans limit and 
optically thin material), the observed flux is 



oc K{v)MdTd, 



(1.2) 



where k(z^) and is the opacity and the temperature of the dust, is the mass of 
the dust in the disk. If can be obtained otherwise, and n{v) is known, the mass of 
the dust disk can be calculated. If the opacity has a power-law dependence {k{v) oc z^^), 
then the flux is Fy oc z^" with a = 2 + (3. Using multi-wavelength measurements, we 
can obtain /3. The /3 parameter in the interstellar rnatter is /3ism = 1.8 it 0.2, but mm- 



observations of disks (most recently bv lRicci et al.l (j2010l )) show that /3 is smaller than 
the ISM value, it is around 0.5- 1. There can be several reasons why the /? parameter is 



reduced in disks ( Drains 



20061 ). Some of the emission might come from optically thick 



regions, therefore the assumption used in deriving Eq. 11.21 is not valid. The chemical 
composition of dust in the disks can be very differen t from that of the ISM dust. Grain 
growth can change the size distribution of particles ([Birnstiel et al.l . |2010| ) . The opacity 
model might not be correct, e.g., the opacity of fractal structures can be quite different 
from non-fractal, spherical particles. 

Millimeter, sub-mm observations, in a similar w ay as IR observations , can b e used to 
infer disk lifetimes. Using these measurements, [Andrews Williams! (j2005l ) obtained 
the same disk lifetime as from IR measurements. 

UV exces s and magnetosphe ric emission lines can be used to determine accretion rates 



of disks (Calvet et al. 



2000). The typical accretion rate for young stars is 10 ^ Mg 







yr ^; for stars of age 1 Myr, it is 10 ^ M© yr ^. As the accretion rate is decreasing in 
time, one caii infer the disk lifetime based on UV observations to be also a few 10^ yrs 
jcalvet et al.l . mm . 



These results of disk lifetime are compelling. The IR excess measures the 'survival 
time' of the small dust grains around 1 AU, the sub-mm measurements are based on °° 

O 

the presence of dust at large radii, while the accretion signature means that gas can be - 
transported onto the surface of the star directly from the gas disk. These three methods 
estimate the disk lifetime measuring entirely different disk material, still they obtain 
the same characteristic timescale. These observations imply that the disk dispersion 
happens across a wide range of radii in a relatively short time. 
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Line emissions of molecular species provide a unique way to determine the physical 
conditions in disks. Using CO observations, the radial and vertical temperature profile of 
disks can be determined. The observation of different CO isotopologues revealed that the 
outer disk is smaller for ^^CO and C^^O than ^^CO, which suggest that photodissociation 
is present at the outer disk. The condensation of CO onto grain surfaces can also 
be observed. These results and the prospects of t he At acama Large Millimeter Array 
(ALMA) are summarized in Guilloteau Sz Dutrey ( 20081 ). 



1.3 Dust experiments in the laboratory 



An ext ensive overview of the available laboratory experiments is given in Chapter 
and in iBlum &: WurmI ([200a) • Here we shortly review the general properties of these 



experiments concentrating on the microphysics of aggregates. 

Monomers are solid bodies which serve as building blocks of aggregates. Often these 
monomers are represented as spheres for simplicity, but in general they can have any 
shape, e.g., ellipsoids or irregular structures. If we assume that the monomers are 
electrically neutral and non-magnetic, these monomers stick together via surfaces forces 
(dielectric van der Walls forces for silicates and molecular dipole interactions for ices). 
The strength of the bond between the n ionomers can be c haracterized by the contact 
force, which can be measured in the lab (jHeim et al.l . Il999l ) and it is given by 



(1-3) 



where 7^ is the specific surface energy and R is the local radius of surface curvature (for 
monomers of different size it is i? = aia2(ai + 02), where ai and 02 are the radii of the 



Poppe et al.l (I2OOOI ) measured the maximum velocity for sticking between 



monomers ) . 

monomers of different sizes impacting onto a smooth surface with different velocities 
and found that this threshold velocity is around 1 m/s for micron sized silicates, and 
that it is decreasing with increasing size. 

At intermediate relative velocities, the collision energy of the aggregates can be dis- 
sipated via restructuring and the two aggregates stick to gether. The most ef f ective 
channel for restructuring is the rolling of two mo n omers (iDoininik &: Tielens 



Wada et al. 



2007 



2008 



Paszun &: Dominik 



20091). 



Heim et al 



1997 



(I1999I ) also measured 



the rolling energy between two monomers, which is the energy needed to roll a monomer 
by 90 degrees, and found that it linearly depends on the local surface curvature of the 
monomers. 

If we want to break a contact by pulling the monomers apart, we have to pull the two 
monomers with a force that is higher than the contact force. The interesting feature of 
the contact area is that it can be stretched further apart such that the distance of the 
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center of mass coordinates of the monomers can actuahy be larger than the sum of the 
radii. The elastic energy stored in this 'neck' is then transformed into elastic waves and 
slowly dissipates. Two mono mers can also twist and sli de, but these motions seem to 
be less important in collisions (jPominik fc Tielend . 119971 ). 



Most of these results were obtained using silicate spheres that are smaller than 1 mi- 
cron in size, thus a simple shape and a mono-disperse size distribution were used. The 
picture is more complicated if we assume a size distribution of monomers f Pominik &: 
Tielens. 119971 ) or irregular monomer shapes (iBluml . l2004l ). Although we have a qualita- 



tive picture how these effects change the particle properties, no exhaustive and general 
investigations were made. There is also uncertainty regarding the properties of dif- 
ferent monomer-materials. It is experimentally not st udied how mi c romet er-sized ice 

([2002) performed 



Kouchi et al 



monomers or monomers with organic mantels behave, 
experiments with millimeter sized grains with organic mantel and found that the sticking 
threshold velocity was several orders of magnitude higher than for same sized silicate 

so increases the sti cking threshold velocity but to a lesser extent 



grains. Frost mantel a 
than organic mantels (IHatzes et al.l . Il99ll ) 



1.4 The importance of dust - theoretical considerations 



The growth of the particles is governed by the Smoluchowski equation, which is intro- 
duced and described in Chapter 12.41 As discussed earlier in this chapter, the particles 
collide and grow because they have a relative velocity in the gas disk. The strength of 
the relative velocity depends on the aerodynamical properties (the stopping time) of the 
aggregates. These properties, as well as the different sources of relative velocities are 
discussed in Chapters 14.2.21 and 15.2.31 

It is clear from the previous section that the observations of disks are strongly influenced 
by the properties of the dust. In this section we review what role the dust plays in various 
aspects of theoretical modeling of disks and planet formation. 

Dust is the main source of opacity where the temperature is below the evaporation 
temperature of the dust (below 1500 K). For an individual dust particle the cross- 
section to radiation at a given wavelength will depend on the particle size, structure 
(fractal or non-fractal, spherical or more complex shape) and composit ion. The total 



opacit y at a given location also depends on the particle size distribution. ISemenov et al 



(|2003l ) describes how to calculate the Rosseland mean opacities for particles as a function ° 
of temperature, assu ming that the particles follow a modified MRN size distribution ° 



(Pollack et al 



19851 ) and that these particles are homogeneous and spheric al. It is also 



possib l e to calculate opa cities of arbitrarily complex aggregates (see e.g. 



Shen et al 



(|2008h 



Min et al. 



(120071 )). However, usually the uncertainty coming from the amount 
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and the size of the dust particles is higher than the uncertainty in opacities. Therefore, 
in most cases approximate opacity formulae are useful. 

The dust also effects the ionization state of the protoplanetary disk. The sources of 
ionization can be thermal (from the central star), and non-thermal (X-rays, cosmic rays, 
^^Al decay). Thermal ionization is effective at the inner disk (< 0.1 AU from the central 
star). X-rays and cosmic rays can ionize a layer of thickness with 100 g/cm? column 
density on both sides of the disk. The decay of ^^Al provides an ionization source that 
is present at every location in the disk. Dust affects ionization in two ways. It acts as 
charg e carrier and it provides surface where electrons and ions can recombine (jOkuzumi , 
200g). 



Although the ionization fraction is almost negligible, it has a very important affect 
for angular momei itum redistribution in disk s. To simulate viscous spreading, usually 
the description by IShakura fc SunyaevI (119731 ) is used, where the turbulent viscosity is 
parameterized. Although this model is widely used, it does not explain what the physical 
source of the turbulent viscosity are. The best candidate-mechanism is currently the 
ma gneto-rotational in s tabili ty (MRI), which was first discussed in the context of disks 
by iBalbus &: Hawleyi (|l99ll ). The basis of the instability is that even low ionization 
fractions make it possible to couple the magnetic field and the gas dynamically. If there 
are two fluid parcels orbiting at different radii, the magnetic field acts as a spring between 
these two fluid parcels. This force causes the inner one to be slowed down by the outer 
one, and the outer one to be sped up by the inner one. As the inner one loses energy 
its orbit shrinks, while the outer one gains energy and moves out. During this process 
the distance, thus the force between these parcels increases further and the process runs 
away. 



here can be a region in the disk (typically between 0.2 and 4 AU - iD'Alessio et al 



(|l998l )). where the i onization , fracti on is smaller than the minimum value required for 



the MRI. Therefore, 



Gammid (jl996l ) proposed a layered accretion disk model, in which a 



'dead zone' is present at the midplane of the disk. The dead zone has suitable properties 
for planet formation. A pressure bump is present at the edges of the dead z ones, where 



dust particles can be accumulat ed and planetesimals can be formed (e.g. iLvra et al 



(|2008h 



Dzvurkevich et al.l (j2O10l ). 



Shear box simulations of MRI-driven turbulence showed that dust particles of decimeter 
in size can be efficiently concentrated in the turbulent eddies. This way planetesimals 
can be formed d i rectly from decimeter sized bodies avoiding all size regimes in between 
(jJohansen et al.l . boOTI ). 



We restricted the discussion on MRI as a possible driving mechanism for angular mo- 
mentum redistributio n within disks, but there are several other can didates, like the 



baroclinic instability (iKlahr fc Bodenheimed . l2003l ). shear instability (iDubrulle et al 
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20051 ) , gr avitational spiral waves (jPickett et al.l . |2003| ) , or global magnetic fields threading 
the disk Jstehle fc Soruitl . boOlh . 



Dust particles are important for chemistry in disks. From the chemical point of view, 
disks can be divided into three regions according to the height above the midplane. 
These a r e the photon dominated layer, warm molecular layer, and the midplane layer 
(|Bergin . 2009 ). In these layers different type of chemistry are dominating, which are 
the gas phase chemistry (at the hot and highly ionized photon dominated layer), gas- 
grain chemistry and grain surface chemistry (i n the colder reg i ons o f the disk). We also 
have to distinguish the inner and outer disks (ISemenov et al.l . l2010l ) . In the inner disk, 
chemical equilibrium is reached within ~ 100 yrs. However, the material in the outer 
disk can have much longer reaction timescales (~ 10^ yrs) due to the low densities and 
temperatures, thus kinetic chemistry must be used to obtain the abundance of different 
chemical species. 



1.5 The outline of the thesis 

This thesis investigates the initial growth of dust particles using a Monte Carlo (MC) 
model which includes a collision model that is based on the currently known laboratory 
experiments on dust aggregates. Chapter [2] describes the basic numerical method which 
is used in this thesis. We test the method by comparing the results against analytical 
solutions and highlight the strengths and weaknesses of the MC method. 

Chapter [3] describes the laboratory based collision model. The assumptions, fitting 
formulas and parameters, and extrapolations that enter the model are discussed in detail. 

Chapter U] combines the results of the previous chapters by incorporating the collision 
model into the MC code. In this chapter, local box simulations are performed at the 
midplane of three different disk models at 1 AU distance from a solar mass star. In this 
chapter we show the importance of bouncing collisions and how dust evolution is altered 
compared to previous simpler collision models. 

Chapter [5] is the extension of the previous OD models. In this chapter we investigate the 
settling and growth of dust particles in a ID vertical column of disks. 

Finally, we review the future prospects of this field in Chapter [H 
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Chapter 2 



A representative particle 
approach to the coagulation of 
dust particles 



Based on 'A representative particle approach to coagulation and fragmentation of dust 
aggregates and fluid droplets' by A. Zsom &: C. P. DuUemond published in A&A, 489, 
931. 



2.1 Introduction 



Dust particle aggregation is a very common process in various astrophysical settings. 
In protoplanetary disks the aggregation of dus t particles forms the very initial step of 
planet formation (see e.g. Dominik et al. 120071 ). It also modifies the optical properties 
of the disk, and it has influence on the c hemistry and free electr on abundance in a disk 
(Sano et al. l2000l : Semenov et al. l2004l : Ilgner & Nelson l200a ). The appearance and 
evolution of a protoplanetary disk is therefore critically affected by the dust aggregation 
process. In sub-stellar and planetary atmospheres the aggregation of dust particles and 
the coagulation of fluid droplets can affect the structure of cloud layers. It can therefore 
strongly affect the spectrum of these objects and influence the local conditions within 
these atmospheres. The process of aggregation/coagulation and the reverse process of 
fragmentation or cratering are therefore important processes to understand, but at the 
same time they are extremely complex. 



Traditional methods solve the Smoluchowski equation for the particle mass distribution 
function f{m), where f{m) is defined such that f(m)dm denotes the number of particles 
per cubic centimeter with masses in the interval [m, m + dm]. This kind of metho d has 



been used i n many pap ers on dust co agulat ion before (e.g. N akaga wa et al 
denschilling 



1981 



1984 



1997: Schmitt et al. 



19971 : Suttner & Yorke 



; Wei- 



199S; Tanaka et al. 



2005 
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Dullemond & Dominik l2005l : Nomura & Nakagawa 120061 ) . Methods of this kind are 



efficient, but have many known problems. First of all a coarse sampling of the particl e 



mass leads to systematic errors such as the acceleration of growth (Ohtsuki et al. Il990l ). 
High resolution is therefore required, which may make certain problems computationally 
expensive. Moreover, if one wishes to include additional properties of a particle, such as 
porosity, charge, composition etc, then each of these properties adds another dimension 
to the problem. If each of these dimensions is sampled properly, this can quickly make 
the problem prohibitively computationally expensive. Finally, the traditional methods 
are less well suited for modeling stochastic behavior of particles unless this stochastic 
behavior can be treated in an averaged way. For instance, in protoplanetary disks if the 
stopping time of a particle is roughly equal to the turbulent eddy turn-over time, then 
the velocity of a particle with respect to the gas is stochastic: at the same location there 
can exist two particles with identical properties but which happen to have different ve- 
locities because th ey ent ered the eddy from different directions (see e.g. the simulations 



by Johansen et al. l2006l ). 



To circumvent problems of this kind Ormel et al. (|2007l ) have presented a Monte Carlo 
approach to coagulation. In this approach the particles are treated as computational 
particles in a volume which is representative of a much larger volume. The simulation 
follows the life of N particles as they collide and stick or fragment. The collision rates 
among these particles are computed, and by use of random numbers it is then determined 
which particle collides with which. The outcome of the collision is then determined de- 
pending on the properties of the two colliding particles and their relative impact velocity. 
This method, under ideal conditions, provides the true simulation of the process, except 
that random numbers are used in combination with collision rates to determine the next 
collision event. This method has many advantages over the tradiational methods. It is 
nearly trivial to add any number of particle properties to each particle. There is less 
worry of systematic errors because it is so close to a true simulation of the system, and it 
is easy to implement. A disadvantage is that upon coagulation the number of computa- 
tional particles goes down as the particles coagulate. Ormel et al. solve this problem by 
enlarging the volume of the simulation and hence add new particles, but this means that 
the method is not very well suited for modeling coagulation within a spatially resolved 
setting such as a hydrodynamic simulation or a model of a protoplanetary disk. 

It is the purpose of this chapter to present an alternative Monte Carlo method which can 
quite naturally deal with extremely large numbers of particles, which keeps the number 
of computational particles constant throughout the simulation and which can be used 
in spatially resolved models. 
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2.2 The method 

2.2.1 Fundamentals of the method 

The fundamental principle underlying the method we present here is to follow the be- 
havior of a limited number of representative particles whose behavior is assumed to be a 
good representation of all particles. In this approach the number of physical particles N 
can be arbitrarily large. In fact it should be very much larger than the number of rep- 
resentative particles n, so that the chance that one representative particle collides with 
another representative particle is negligible compared to the collisions between a repre- 
sentative particle and a non-representative particle. In other words, if ^ n, we only 
need to consider collisions between a representative particle and a non-representative 
particle. The number of collisions among representative particles is too small to be sig- 
nificant, and the collisions among non-representative particles are not considered because 
we focus only on the behavior of the representative particles. 

Suppose we have a cloud of dust with N = 10^^ physical pa rticles , with a specific size 



distribution, for instance, MRN (Mathis, Rumpl & Nordsieck 1 1 9771 ) . Let the total mass 
of all these particles together be Mtot and the volume be V. We randomly pick n 
particles out of this pool, where n is a number that can be handled by a computer, for 
instance, n = 1000. Each representative particle i has its own mass rrii and possibly 
other properties such as porosity pi or charge Cj assigned to it. We now follow the life of 
each of these n = 1000 particles. To know if representative particle i = 20 collides with 
some other object, we need to know the distribution function of all physical particles 
with which it can collide. However, in the computer we only have information about 
the n representative particles. We therefore have to make the assumption that the 
distribution function set up by the n representative particles is representative of that of 
the A'^ physical particles. We therefore assume that there exist 

nrukV 

physical particles per cubic cm with mass mk, porosity pk, charge Ck etc., and the same 
for each value of k, including k = i. In this way, by assumption, we know the distribution 
of the N physical particles from our limited set of n representative particles. One could 
say that each representative particle represents a swarm of Mtot/nrrii physical particles 
with identical properties as the representative one. One could also say that the true 
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distribution of particles is, by assumption, that of the n representative ones. The % b b b o 
rate of collisions that representative particle i has with a physical particle with mass 5 2 2 5 5 
etc. is then: 

ri{k) = UkaikAvik = ^(JikAvik (2.2) 
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where is the cross-section for the collision between particles with properties i and '^'^ j^^B '° g 
k, and Avik is the average relative velocity between these particles. The total rate of □ ] 2 
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collisions that representative particle i has with any particle is then: 

n = J]ri(fc) (2.3) 

k 

and the total rate of collisions of any representative particle is 

r = J2ri (2.4) 



The time-evolution of the system is now done as follows. Let to be the current time. We 
now randomly choose a time step St according to: 

St = — - log(ran(seed)) (2.5) 

where ran(sccd) is a random number uniformly distributed between and 1. This means 
that a collision event happens with one of the representative particles at time t = to + 6t. 
The chance P{i) that the event happens to representative particle i is: 

P{{) = ^ (2.6) 

So we can choose, using again a random number, which representative particle i has 
undergone the collision event. We now need to determine with what kind of physical 
particle it has collided. Since the distribution of physical particles mirrors that of the 
representative ones, we can write that the chance this particle has collided with a physical 
particle with properties k is: 

ptkli) = ^ (2.7) 

With another random number wc can thus determine which k is involved in the collision. 
Note that k can be i as well, i.e. the representative particle can collide with a physical 
particle with the same properties, or in other words: a representative particle can collide 
with a particle of its own swarm of physical particles. 

Now that we know what kind of collision has happened, we need to determine the out- 
come of the collision. The most fundamental part of our algorithm is the fact that only 
representative particle i will change its properties in this collision. Physical particle k 
would in principle also do so (or in fact becomes part of the new representative particle), 
but since we do not follow the evolution of the physical particles, the collision will only 
modify the properties of representative particle i. By assumption this will then automat- 
ically also change the properties of all physical particles associated with representative 
particle i. Statistically, the fact that the particles k are not modified is "corrected for" 
by the fact that at some point later the representative particle k will have a collision 
with physical particle i, in which case the properties of the k particles will be modified 
and not those of i. This then (at least in a statistical sense) restores the "symmetry" 
of the interactions between i and k. If the collision leads to sticking, then the resulting 
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particle will have mass m = mi + rrik. This means that representative particle i will from 
now on have mass rrii ■i^ rrii + m^. Representative particle k is left unaffected as it is 
not involved in the collision. The interesting thing is now that, because by assumption 
the representative particle distribution mirrors the real particle distribution, the swarm 
of physical particles belonging to the modified representative particle i now contains 
fewer physical particles, because the total dust mass M = Mt^t/n of the swarm remains 
constant. 

If a collision results in particle fragmentation, then the outcome of the collision is a 
distribution function of debris particles. This distribution function can be written as a 
function /d(m) of debris particle mass, such that 

mfd{m)dm = nii + rrik (2-8) 







and the function fd{fn) has to be determined by laboratory experiments or de tailed 



2007 for a 



computer simulations of individual particle collisions (see Dominik et al. 
review). The new value of rrii for the representative particle is now randomly chosen 
according to this distribution function by solving the equation 

j-ifi 

/ mfd{m)dm = ran(seed)(mi + ruk) (2.9) 
Jo 

for rh and assigning ^ fh. In other words: we randomly choose a particle mass 
from the debris mass distribution function, i.e. the choice is weighed by fragment mass, 
not by fragment particle number. This can be understood by assuming that the true 
representative particle before the collision is in fact just a monomer inside a larger 
aggregate. When this aggregate breaks apart into for instance one big and one small 
fragment it is more likely that this representative monomer resides in the bigger chunk 
than in the smaller one. 

After a fragmenting collision the mi will generally be smaller than before the collision. 
This means that the number of physical particles belonging to representative particle i 
increases accordingly. Note that although the collision has perhaps produced millions 
of debris particles out of two colliding objects, our method only picks one of these 
debris particles as the new representative particle and forgets all the rest. Clearly if 
only one such destructive collision happens, the representative particle is not a good 
representation of this entire cloud of debris products. But if hundreds such collisions 
happen, and are treated in the way described here, then the statistical nature of Eq. ()2.9p 
ensures that the debris products are well represented by the representative particles. 

The relative velocity Av can be taken to be the average relative velocity in case of random 
motions, or a systematic relative velocity in case of systematic drift. For instance, for 
Brownian motion there will be an average relative velocity depending on the masses 
of both particles, but differential sedimentation in a protoplanetary disk or planetary 
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atmosphere generates a systematic relative velocity. Also, for the Brownian motion or 
turbulent relative velocity one can, instead of using an average relative velocity, choose 
randomly from the full distribution of possible relative velocities if this is known. This 
would allow a consistent treatment of fluctuations of the relative velociti es wh ich could 
under some circumstances become important (see e.g. Kostinski & Shaw [20051). 



2.2.2 Computer implementation of the method 

We implemented this method in the following way. For each of the representative parti- 
cles we store the mass rrii and all other properties such as porosity, charge, composition 
etc. Before the start of the Monte Carlo procedure we compute the full collision rate 
matrix ri{k), and we compute the rj as well as r. For these collision pairs {i, k) we now 
have to determine the cross section of particles as well as their systematic relative ve- 
locity, such as different drift speeds, and the random relative velocity, such as Brownian 
or turbulent motion. The random motions can be determined with a random number 
from the relative velocity probability distribution function if that is known. If that is 
not known in sufficient detail, one can also take it to be the average relative velocity, for 
which more often analytic formulae exist in the literature. 

We determine beforehand at which times igav.n we want to write the resulting mj and 
other parameters to a file. The simulation is now done in a subroutine with a do- 
while loop. We then determine 5t using a random number (see Eq. 12. Sp . and check if 
t + 6t < tsav.n, where tsav,n IS the next time when the results will have to be stored. If 
t + 6t < isav.nj then a collision event occured before tsav,n- We will handle this event 
according to a procedure described below, we set t <— t + 6t and then return to the 
point where a new 6t is randomly determined. If, on the other hand, t + 6t > tsav,n 
then we stop the procedure, return to the main program and set t ^ isav.n- The main 
program can then write data to file and re-call the subroutine to a time tsav,n+i or stop 
the simulation altogether. Note that when the subroutine is called again for a next time 
interval, it does not need to know the time of the previously randomly determined event 
which exceeded isav,n- Of course, one could memorize this time and take that time as 
the time of the next event in the next time interval, but since the events follow a Poisson 
distribution, we do not need to know what happened before tsav,n to randomly determine 
the new time t + St of the next event. 

Now let us turn to what happens if a collision event occurs, i.e. occurs between time t 
and t+6t. We then first determine which representative particle i is hit, which is done by 
generating a random number and choosing from the probability distribution of collision 
rates, as described in Section [2.2.11 Similarly we determine the non-representative par- 
ticle with which it collides, or in other words: we determine the index k of the "swarm" 
in which this non-representative particle resides. Finally, we must determine the impact 
parameter of the collision, or assume some average impact parameter. 
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Now we employ a model for the outcome of the collision. This is the collision subroutine 
of our Monte Carlo method. It is here where the results of laboratory experiments 
come in, and the translation of such experiments into a coagulation kernel is a major 
challenge which we do not cover here. The collision model must be a quick formula 
or subroutine that roughly represents the outcomes of the detailed laboratory collision 
experiments or detailed numerical collision models. It will give a probability function 
fd for the outcoming particle masses and properties. Prom this distribution function 
wc pick one particle, and from this point on our representative particle i will attain 
this mass and these properties. The collision partner k will not change, because it is a 
non-represetative particle from that swarm that was involved in the collision, and wc do 
not follow the life of the non-representative particles. We therefore ignore any changes 
to that particle. 

We now must update rj{l) for all / with fixed j = i and for all j with fixed I = k: we 
update a row and a column in the rj(l) matrix. Having done this, we must also update 
rj for all j. This would be an process, which is slow. But in updating rj{l) we know 
the difference between the previous and the new value, and we can simply add this 
difference to rj for each j. Only for j = i we must recompute the full rj again, because 
there all elements of that row have been modified. Using this procedure we assure that 
we limit the computational effort to only the required updates. 



2.2.3 Acceleration of the algorithm for wide size distributions 

One of the main drawbacks of the basic algorithm described above is that it can be very 
slow for wide size distributions. Consider a swarm of micron sized dust grains that are 
motionless and hence do not coagulate among each other. Then a swarm of meter sized 
boulders moves through the dust swarm at a given speed, sweeping up the dust. Let us 
assume that also the boulders are not colliding among each other. The only mode of 
growth is the meter-sized boulders sweeping up the micron sized dust. For the boulder 
to grow a factor of 2 in mass it will have to sweep up 10^^ micron sized dust particles. 
Each impact is important for the growth of the boulder, but one needs 10^® such hits to 
grow the boulder a factor of 2 in mass. The problem with the basic algorithm described 
above is that it is forced to explicitly model each one of these 10^® impacts. This is 
obviously prohibitively expensive. 

The solution to this problem lies in grouping collisions into one. Each impact of a dust , 
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grain on a boulder only increases the boulder mass by a minuscule fraction. For the % 'o b L o 
growth of the boulder it would also be fine to lower the chance of an impact by 10^®, S 2 S 5 S 



but if it happens, then 10 particles impact onto the boulder at once. Statistically 
this should give the same growth curve, and it accelerates the method by a huge factor. 
However, it introduces a fine-tuning parameter. We must specify the minimum increase 
of mass for coagulation {drumax)- If we set drumax to, for instance, 10%, then we may 
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expect that the outcome also has errors of the order of 10%. This error arises because 
by increasing the mass of the bigger body in steps of 10%, we ignore the fact that the 
mass at some time should in fact be somewhere in between, which cannot be resolved 
with this method. This is, however, not a cumulative error. While the mass of the 
bigger body may sometimes be too low compared to the real one, it equally probably 
can be too large. On average, by the Poisson nature of the collision events, this averages 
out. But it is clear that the smaller we take this number, the more accurate it becomes 
- but also the slower the method becomes. It is therefore always a delicate matter to 
choose this parameter, but for problems with a large width of the size distribution this 
acceleration is of vital importance for the usability of the method. 



2.2.4 Including additional particle properties 

We mentioned briefly the possibility of adding more particle properties to each represen- 
tative particle. This is very easy to do, and it is one of the main advantages of a Monte 
Carlo method over methods that directly solve the integral equations for coagulation. 
One of the main properties of interest to planet formation is porosity or fractal structure 
of the aggregate. Two aggregates with the same mass can have vastly different behavior 
upon a collision if they have different compactness. A fluffy aggregate may break apart 
already at low impact velocities while a compact aggregate m ay sim ply bounce. Upon 



collisions these properties may in fact also change. Ormel et al. (|2007l ) studied the effect 
of porosity and how it changes over time, and they also used a Monte Carlo approach 
for it. 

If one wishes to include particle properties in a traditional method which solves the 
integral equations of coagulation (the Smoluchowski equation), then one increases the 
dimensionality of the problem by 1 for each property one adds. With only particle mass 
one has a distribution function f(m,t) while adding two particle properties pi and p2 
means we get a distribution function f{m,pi,p2,t), making it a 4-dimensional problem. 
Methods of this kind must treat the complete phase space spanned by {m,pi,p2,t). This 
is of course possible, but computationally it is a very challenging task (see Ossenkopf 



19931 ). In contrast, a Monte Carlo method only sparsely samples phase space, and it 
samples it only there where a significant portion of the total dust mass is. A Monte 
Carlo method focuses its computational effort automatically there where the action is. 
The drawback is that if one is interested in knowing the distribution function there 
where only very little mass resides, then the method is inaccurate. For instance, in a 
protoplanetary disk it could very well be that most of the dust mass is locked up in big 
bodies (larger than 1 meter) which are not observable, and only a promille of the dust 
is in small grains, but these small grains determine the infrared appearance of the disk 
because they have most of the solid surface area and hence most of the opacity. In such 
a case a Monte Carlo method, by focusing on where most of the mass is, will have a very 
bad statistics for those dust grains that determine the appearance of the disk. For such 
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Figure 2.1: Results of the test problem with N swarms of small particles and N 
swarms of big bodies, as discussed in Section USUI Left: histogram of the final masses 
of the bodies relative to the initial mass of the big bodies, for N = 500. Right: average 
mass relative to the initial mass of the big bodies as a function of N. 



goals it is better to use the traditional methods. But if we are interested in following 
the evolution of the dominant portion of the dust, then Monte Carlo methods naturally 
focus on the interesting parts of phase space. 



2.3 Discussion of the method 



2.3.1 Conservation of particle number 



There are a few peculiarities of the method described here that may, at first sight, appear 
inconsistent, but are statistically correct. For instance if we return to our example of a 
swarm of tiny particles and a swarm of boulders, i.e. n = 2 with representative particle 1 
being a micron sized particle and representative particle 2 being a meter sized particle, 
then we encounter an apparent paradox. We again assume that collisions only take place 
between 1 and 2, but not between 1 and 1 or 2 and 2. The chance that representative 
particle 1 hits a meter size particle is much smaller than the chance that representative 
particle 2 hits a micron size particle. What will happen is that representative particle 
2 will have very many collision events with small micron size grains, and thereby slowly 
and gradually grows bigger, while representative particle 1 will only have a collision 
with big particle after a quite long time and immediately jumps to that big size. While 
representative particle 2 grows in mass, the number of big physical particles decreases 
in order to conserve mass. This may seem wrong, because in reality the number of big '^tisusp ssom mjou 

oil II 

boulders stays constant, and these boulders simply grow by sweeping up the small dust. ° ° ° ° 

p IT) ^ c 

The solution to this paradox is that the average time before representative particle 1 hits 
a big {k = 2) particle is of the same order as the time it takes for representative particle 
2 to grow to twice its mass by collecting small particles. So, very roughly, by the time 
the big particle has doubled its mass, and therefore the number of physical particles 
belonging to k = 2 has reduced by 50%, the representative particle 1 has turned into 
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a big particle, corresponding, statistically, to the other 50% of big particles that was 
missing. If we are a bit more precise, the statistics do not add up precisely in this way 
if we have only 1 swarm of small and 1 swarm of big bodies. If, however, one has N 
swarms of small and N swarms of big bodies, and again assume that only the big bodies 
can sweep up the smaller ones, then if ^ 1 the statistics adds up perfectly: one finds 
that after all the growth has taken place, the average mass of the bodies is twice that 
of the original big bodies. In Figure 12.11 we do precisely this experiment, and the left 
panel shows that for N = 500 the mass distribution of the big bodies averages to the 
right value, albeit with a spread of 10% FWHM while in reality this spread should be 
0. The right panel shows how the average final mass depends on A'^. For small the 
statistics clearly do not add up, but for large they do and produce the right value 
(final mass is twice initial mass of the big bodies). So statistically the number of big 
particles is restored to the correct value, but there is then unfortunately still a large 
statistical noise on it. The particle number is therefore not exactly conserved in our 
method, but statistically it is. 



2.3.2 The number of representative particles 

It is obvious that for high number of representative particles n we will get better results 
than for low n. But there are two issues here. First of all, the higher n, the better 
the representative particles represent the true physical distribution of particles. For 
problems that result in wide size distributions this is all the more crucial. An inaccurate 
representation of the true size distribution could lead to systematic errors. But another 
reason for taking a high n is simply because we want our end-result to have as little as 
possible noise. If the result is too noisy, then it is useless. Taking n too big, however, 
makes the code slow because more representative particles have to be followed, and for 
each of these particles we must check for a larger number of possible collision partners k. 
The problem scales therefore as v?. If the expected size distributions are not too wide, 
one can use an intermediately large number of representative particles, say n, for the 
simulation, but redo the simulation m times such that n = mn, a nd ave rage the results 



of all m simulations. This approach was also used by Ormel et al. (120071 ). This gives the 
same amount of noise on the end-result, but scales as n'^m = n'^/m, which is m times 
faster than the n? scaling. This works, however, only if the coagulation/fragmentation 
kernel is not too sensitive to the exact distribution of collision partners. 

Interestingly, if the kernel is very insensitive to the exact distribution of collision part- 
ners, then, in principle, one could run the model with only a single representative particle 
n = 1, because the collision partner of representative particle i could be equal to k = i. 
Of course, a single representative particle means that we assume that all physical par- 
ticles have the same size, or in other words: that we have an infinitely narrow size 
distribution. 
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To decide about the sufficient number of representative particles, one has to compare the 
results of the MC code with the analytical solutions of the three test kernels (see Section 
12. 4p . In a given time of the simulation the mean mass and the shape of the distribution 
function for all three test kernels must be followed accurately. It is especially important 
to reproduce the linear and product kernels accurately as the realistic kernels of dust 
particles are similar to these. 

Of course the progression from 'not sufficient' and 'sufficient' number of representative 
particles is smooth and in general the more representative particles we use, the more 
accurate the produced result will be. The sufficient number of representative particles 
(n) as given in Section [2.4p are only suggestions, the error of the distribution functions 
were not quantified. 



2.3.3 Limitations of the method 



One of the fundamental limitations of the method described here is that we assume 
N ^ n. We can model the growth of particles by coagulation in a protoplanetary 
disk or in a cloud in a planetary atmosphere, but we can not follow the growth to the 
point where individual large bodies start to dominate their surroundings. For instance, 
if we wish to follow the growth of dust in a protoplanetary disk all the way to small 
planets, then the method breaks down, because A'^ is then no longer much bigger than 
n, and interactions among representative particles become likely. Also, for t he same 
reason, run-away growth problems such as electrostatic gelation (Mokler 120071 ) cannot 
be modeled with this method. 



Another limitation is encountered when modeling problems with strong growth and 
fragmentation happening at the same time. This leads to very wide size distributions, 
and the typical interval between events is then dominated by the smallest particles, 
whereas we may be interested primarily in the growth of the biggest particles. In such a 
situation a semi-steady-state can be reached in which particles coagulate and fragment 
thousands of times over the life time of a disk. The Monte Carlo method has to follow 
each of these thousands of cycles of growth and destruction, which makes the problem 
very "stiff". Methods using the integral form of the equations, i.e. the Smoluchowski 
equation, can be programmed using implicit integration in time so that time steps can 
be taken which are much larger than the typi cal tim e scale of one growth-fragmentation 
cycle without loss of accuracy (Brauer et al. |2008| ). This is not possible with a Monte 
Carlo method. 
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Figure 2.2: Test against the constant kernel (i^i.j — 1). The particles were 
binned and the distribution function was produced at dimensionless times t = 
0, 10°, 10\ 10^, 10^, 10'*, 10^. The dashed lines show the analytical solution. This run 
was produced by simulating 200 representative particles five times and producing the 
average of these. In this case dnimax is 0.1. 



2.4 Standard tests and results 



In this section we test our coagulation model with kernels that have analytical solutions. 
Furthermore we show the first results of applying this model to protoplanetary disks 
introducing Brownian motion and turbulence induced relative velocities as well as a new 
prop erty of dust particles namely the porosity (or enlargement factor, see Ormel et al 



20071 ). and a simple fragmentation model. 



To follow dust coagulation and fragmentation, one has to follow the time evolution of the 
particle distribution function at a given location in the disk (/(y,t)), where y contains 
the modeled properties of the dust grains, in our case these will be the mass (m) and 
the enlargement factor {^), f{m,^,t). 

In most of the coagulation models so far the only used dust- property was t h e par ticle 



mass. Then one can use the so called Smoluchowski equation (jSmoluchowskil ()1916l )) to 
describe the time-evolution of f{m): 



df{m) / dm'K{m,m')f{m')+ 



at 



dm'K{m' , m — m')f{m')f{m — m). (2.10) 



The first term on the right hand side represents the loss of dust in the mass bin m 
by coagulation of a particle of mass m with a particle of mass m' . The second term 
represents the gain of dust matter in the mass bin m by coagulation of two grains of 
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Dimensionless moss 

Figure 2.3: Test against the hnear kernel {Kij = mi + nij). The particles were binned 
and the distribution function was produced at dimensionless times t = 0, 4, 8, 12, 16, 20. 
The dashed lines show the analytical solution. This run was also produced by simulating 
200 representative particles five times and producing the average of these. In this case 

drrimax is 0.1. 



mass m' and m — m'. K is the coagulation kernel, it can be written as 



K{mi,m2) = ac{mi,m2) x Af(mi,m2), 



(2.11) 



the product of the the cross-section of two particles and their relative velocity. We 
consider all the three kernels for which there exist analytical solutions: The constant 
kernel [Ki^j = 1), the linear kernel {Kij = rrii + rrij) and the product ker i iel {K ^ 



rrij X m,- 



Wetherih (|l99nl ). 



The analytical solutions are described e.g. in 



Ohtsuki et al. 



19901) and 



We test our method against these three kernels, leaving the enlargement factor un- 
changed, always unity. Further important properties of the dust particles, such as 
material density and volume density, are also always unity. The (dimensionless) time 
evolution of the swarms is followed and at given times the particles are binned by mass 
so that we can produce /(m). On Figures and [23] the y axis shows fjrn ) x m? , the 



mass dens ity pe r bin. The analytical solutions, taken from lOhtsuki et al.l ([19901) and 
Wetherill (jl99ol ). are overplotted with dashed line. The number of particles were chosen 
to be n = 200, m = 5, so altogether 1000 representative particles were used in the model 
except for the product kernel where more representative particles were used to achieve 
better results. 

In the case of the constant kernel (Figure [2. 2p . we started our simulation with MRN size 
distribution (n(a) oc a"^'^), the results were saved at t = 0, 10°, 10\ 10^, 10^, 10^, 10^ 
It is interesting to note that this kernel is not sensitive to the initial size distribution. 
As the system evolves, it forgets the initial conditions. Another interesting property of 
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Figure 2.4: Test against the product kernel {Kij — nii x rrij). The parti- 
cles were binned and the distribution function was produced at dimenslonless times 
t = 0.4,0.7,0.95. The dashed lines show the analytical solution. This run was pro- 
duced by simulating 1000 representative particles ten times and producing the average 
of these. In this case dm„iax is set to be 0.05. 



this kernel that our model can reproduce the analytical solution even with very limited 
number of representative particles (even for n = 5!) but of course with higher noise. It 
is possible to use only one representative particle, which means that the representative 
particle collides with particles from its own swarm which basically results in pure CCA 
growth (Cluster-Cluster Aggregation). Interestingly, the mean mass of the distribution 
function is followed correctly but the shape of the function changes, additional spikes 
appear on it. 

The linear kernel is known to be more problematic because the mean mass of the particles 
grows exponentially with time. Our model, however reproduces this kernel very well, 
too, as it can be seen in Figure 12.31 The results were saved at t = 0, 4, 8, 12, 16, 20. 
We note that using low number of representative particles with this kernel also works 
relatively well, the minimum number of swarms needed to reproduce the exponential 
time evolution of the mean mass is n w 100. This is larger than for the constant 
kernel. It shows that for the linear kernel collisions between particles of unequal mass 
are contributing significantly to the growth, whereas for the constant kernel the growth 
is dominated by collisions between roughly equal size particles. Using n <^ 100 results 
in distorted distribution function: neither the mean mass nor the actual shape of the 
distribution function is correct. 

The product kernel is the hardest to reproduce. The peculiarity of this kernel is the 
following: Using dimenslonless units, a 'run-away' particle is produced around t = 1, 



which collects all the other particles present in the simulation (Wetherill [19901) ■ The 
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difficulty arises in our Monte Carlo code when the mass of the representative 'run- 
away' particle reaches the mass of its swarm. In other words, the number of physical 
particles belonging to the representative 'run-away' particle is close to unity. In this case 
the original assumption of our method (we only need to consider collisions between a 
representative particle and a physical particle) is not valid anymore. However, as Figure 
12.41 shows, we can relatively well reproduce this kernel before t = 1. In the case of this 
kernel, we need approximately n ~ 500 representative particles to correctly reproduce 
it. 

The required CPU time for these test cases is very low, some seconds only. 

We conclude that our Monte Carlo method reproduces the constant and linear test 
kernels without any problem even with low number of representative particles. On the 
other hand the method has difficulties with the product kernel, but before the formation 
of the 'run-away' particle, we can reproduce the kernel. The relatively low number of 
representative particles needed to sufficiently reproduce the test kernels is very important 
for future applications where whole disk simulations will be done and there will likely 
be regions containing low numbers of particles. 



2.5 Applications to protoplanetary disks 

We use the Monte Carlo code to follow the coagulation and fragmentation of dust par- 
ticles in the midplane of a protoplaneta ry disk at 1 AU from the central star. Our disk 
model is identical with the one used by iBrauer et al.l (j2007l ). We proceed step by step. 
First relative velocities induced by Brownian motion and turbulence without the effects 
of porosity are included (Sec. I2.5.ip . 

The next step is to include a fragmentation model (Sec. I2.5.2p . 



In the final step poro sity is included (Sec. I2.5.3p . We use the porosity model described 
in Orr nel et al. (|2007l ). At this point we compare and check again our code with Ormel 
et al. (j2007l ) using their input parameters but not including the rain out of particles. 



2.5.1 Relative velocities 



We include two processes in calculating the relative velocities: Brownian motion and X^ISUSp SSDLU 'LUJOU 
turbulence. ° ° ° ° ° 



Brownian motion strongly depends on the mass of the two colliding particles. The 
smaller their masses are, the more they can be influenced by the random collisions with 
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the gas molecules/atoms. One can calculate an average velocity given by 



Avs{rn^,rn,) = J'J^^^^^^l±^. (2.12) 
y 7rmim2 

For micron sized particles, relative velocity can be in the order of magnitude of 1 cm/s, 
but for cm sized particles this value drops to 10~^ cm/s. If growth is only governed by 
Brownian motion, it leads to very slow coagulation, a narrow size distribution and fluffy 
dust particles, so called cluster-cluster aggregates (CCA). 

The gas in the circumstellar disk is turbulent, thus the dust particles experience ac- 
celeration from eddies with dif ferent sizes and turnover times. This process is very 



complex, but Ormel and Cuzzi (|2007l ) provided limiting closed-form expressions for av- 
erage relative turbulent velocities between two dust particles. Their results are also 
valid for particles with high Stokes numbers. They distinguished three regimes: a.) the 
stopping times of both dust particles are smaller than the smallest eddy-turnover time 
(ti,t2 < t^rj, tightly coupled particles); b.) the stopping time is between the smallest 
and largest turnover time {tri < ti < ti, intermediate regime); c.) the stopping time is 
bigger tha n the largest turnover time (ti > t^, heavy particles). For details see Ormel 



and Cuzzi (120071 ). We used a = 10 for the turbulence parameter. 



To illustrate the relative velocity of dust particles without the effects of porosity, we 
provide Figure 12.51 This contour plot includes Brownian motion and turbulent relative 
velocities. The Brownian motion is negligible for particles bigger than 10"^ cm. 



2.5.2 Fragmentation model 

The collision energy of the particles is 

where fi is the reduced mass. We need to define some quantities of the dust particles. 
Eroii is the rolling ener gy of two monomers. For m onomers of the same size it is given 
by (Dominik & Tielens llM^; Blum k Wurmbood) 



Eroii = ^-^aoFroii, (2.14) 



where qq is the monomer radius, Froii is the rolling force measured by Heim et al. (|l999l ). 
Its value is Fj-oii = (8.5 ±1.6) x 10~^ dyn for Si02 spheres. 

The fragmentation energy is then defined as follows: 

Efrag = Nc X Ebreak ^ 3N X Erolh (2.15) 
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Figure 2.5: The relative velocity caused by Brownian motion and turbulence for dif- 
ferent sized particles. The black line shows the fragmentation barrier. Collision events 
situated between these two lines result in fragmentation if porosity is not included. 
Physical parameters of the disk: the distance from the central star is 1 AU, temper- 
ature is 200 K, the density of the gas is 8.73 x 10^^ particlc/cm'^, and the turbulent 
parameter, a = 10^'^. Parameters of the dust: monomer radius is ag — 0.4/im, material 

density is p = 1.6 g/cm'^. 



where Nc is the total number of contact surfaces between monomers (for simplicity it is 
taken to be 3N, where N is the number of monomers in the particle), Ef,reak is the energy 
needed to break the bond between two monomers (its order of magnitude is similar to 
Eroii for these parameters). 

If the collision energy of two particles is higher than the corresponding fragmentation 
energy, then the aggregate is destroyed and monomers are produced. Note that although 
assuming a complete destruction of the collided dust particles, we are interested in the 
critical energy where the first fragmentation event happens. This is the reason why the 
fragmentation energy is assumed to be lower than the energy needed for catastrophic 
fragmentation. It is a simplification of the model to assume that the debris particles 
will be monomers. Thi s is a very simplified fragmentation model used previously by 



et al. (120071) 



Dull einond &: Dominik (120051). A more realistic model would be the one used by Brauer 



We show the fragmentation barrier in Figure 12.51 with black lines. If collision happens 
in the regime between these two lines, that results in fragmentation. 



2.5.2.1 Results 



A simulation was made including these effects in a specified location of the disk. We 
choose t he location to be 1 A U distance from the central solar type star. Using the disk 
model of . 



Brauer et al 



(j2007l ). the temperature at this distance is approximately 200 K, 
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the density of the gas is 8.73 x 10^^ cm ^, the gas-to-dust ratio is 100 and we choose 
the turbulent parameter to be a = 10~^, the Reynolds number is Re = (based on 



Ormel &: Cuzzi 120071 ) . The dust monomers have the fohowing properties: the monomer 
radius is oq = 0.4 ^um, material density is p = 1.6 g/cm^. With the used parameters 
the fragmentation velocity is Avf^ag ~ 8 m/s, though it is somewhat larger for equal 
sized agglomerates. It is important to note that this value is very sensitive to the 
monomer radius (ao) and material density (p), because smaller/lighter monomers mean 
more contact surfaces (higher N for the same mass) and therefore higher fragmentation 
energy. 

Using these input parameters we simulated the evolution of the dust particles for 3 x 10^ 
years so that we reach an equilibrium between coagulation and fragmentation. Figure [22] 
shows the resulting normalized size distributions in times after t = 3 x W^, 3 x 10^, 
3 X 10^ and 3 x 10'^ years. We used n = 100 particles averaging over m = 100 times (10^ 
particles altogether). The required CPU time to perform this simulation is 1.5 hours 
approximately, dnimax is set to be 0.001 from now on in every simulation. We would 
like to note that giving dnimax (Section 12.2.31) a higher value would decrease the CPU 
time. 

One can see that coagulation happens due to Brownian motion in the beginning of 
the simulation (until 3 x 10^ years) but after that turbulence takes over and the first 
fragmentation event happens after roughly 10^ years. After this event the "recycled" 
monomers start to grow again, but as we see in Figure [231 particles can not reach bigger 
sizes than 0.07 cm. 

We would like to draw attention to the sudden decrease of particles around 0.002 cm in 
Figure [TBI This is the result of the turbulent relative velocity model used here (discussed 
in Sect. 12. 5. ID . At this point the particles leave the 'tightly coupled particles' regime 
and enter the 'intermediate' regime. But the transition in relative velocity between 
these regimes is not smooth, there is a jump in relative velocity from ~20 cm/s to ~60 
cm/s. As a result, particles coagulate suddenly faster and leave this part of the size 
distribution rapidly. Similar 'valleys' can be seen in the following figures with porosity, 
but the feature is less distinct as the stopping times can be different for particles with 
same mass. 



2.5.3 Porosity 



To be able to quantitatively discuss the effect of porosity, we hav e to define the enlarge- 
ment parameter following the discussion of Ormel et al. (|2007l ). If V is the extended 
volume of the grain and V* is the compact volume, than one can define the enlargement 
parameter (^) as 
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Figure 2.6: The evolution of dust particles including the effects of Brownian motion 
and turbulence. Porosity is not included in this model. The particle distribution is 
saved after t = 3 x 10° years - dash-dot line, 3 x lO""^ years - dashed line, 3 x 10^ years 
- dotted line, and 3 x lO*^ years - continuous line. 
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Figure 2.7: The evolution of dust particles including the effects of Brownian motion 
and turbulence. Porosity is included in this model! The x axis shows the compact 
radius. The particle distribution is saved after t = 3x 10° years - dash-dot line, 3 x 10^ 
years - dashed line, 3 x 10^ years - dotted line, and 3 x 10"^ years - continuous line. Note 
that the scaling of the x axis is different from Figure [ 
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V* 



(2.16) 



Compact volume is the volume occupied by the monomers not taking into account 
the free space between the monomer spheres. One can think of it as melting all the 
monomers into a single sphere, the volume of this sphere is the compact volume. We 
use compact radius later on, which is the radius of this sphere. In the previous section 
the mass/volume ratio was constant for the particles. Therefore we could automatically 
calculate the mass of the particle if the radius was known or vice- versa. But from now 
on a particle with given mass m can have a wide range of effective radii depending on 
its enlargement parameter. 

It is essential to know how the enlargement parameter changes upon collisions. We have 
to refine our fragmentation model and intro duce t wo more regimes regarding to collision 
energy. We use the model of Ormel et al. (120071 ) and we only summarize their model 
here. 

The first regime is the low collision energy regime, where the collision energy is smaller 
than the restructuring energy {E < Erestr, where Ej-estr = 5Ej.oii), meaning that the 
particles stick where they meet, the internal structure of the grain does not change. 

The recipe for the resulting enlargement factor after the collision of two particles assum- 
ing that mi > 7712 then is 



7712*1' 5 



add: 



where {^)m is the mass averaged enlargement factor of the colliding particles: 

mi^'i + m2^2 



(2.17) 



nil + ^2 



(2.18) 



Furthermore SccA is the CCA-char acteri stic exponent calculated by detailed numerical 
studies such as Paszun & Dominik (j2006l ) {5ccA = 0.95). ^'a^^jsa necessary additional 
factor for the enlargement factor (for details see Ormel et al. 120071 ): 



'^add = — Wiexp 
mi 



lOmo 



(2.19) 



where m-o is the monomer mass. 



The second regime is the regime of compaction. The internal structure of the monomers 
inside the particle changes, this causes a decreasing porous volume. If the collision 
energy Erestr < -E < Ej^ag: we talk about compaction. In this case the porosity after 
the collision becomes 

* = (l-/c)((^')^-l) + l, (2.20) 
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where fc = E/{NErou) = —AV/V is the relative compaction. One can see that fc 
has to be smaller than unity otherwise ^ in Eq. 12.201 becomes less than unity. But it 
can theoretically happen that E > NEroii- In this long as the total collision 

energy remains below the fragmentation threshold, we assume that after compaction 
this excess energy goes back into the kinetic energy of the two colliding aggregates. 
The two aggregates therefore compactify and bounce, without exchanging mass or being 
destroyed. Bouncing is therefore included in this model, albeit in a crude way. 

The third regime is fragmentation as it was discussed in the previous section {E > Ef^ag)- 
We use the same fragmentation model as before so the result of a fragmenting collision 
are monomers. 



2.5.3.1 Results 



We performed a simulation with exactly the same initial conditions as in the last section 
but we included the porosity as an additional dust property in the model. The result 
can be seen in Figure 12.71 (the required CPU time here is also 1.5 hours). One can 
immediately see that including porosity increases the maximum particle mass by two 
orders of magnitude (five times lar ger p articles in radius). This was already expected 
based on the work of Ormel et al. (j2007l ). although due to rain out of bigger particles, 
they did not simulate particles bigger than 0.1 cm. 

We provide Figure [THl to give an impression how the porosity of the agglomerates change 
during the simulation. The x axis is the compact radius of the particles, the y axis is 
the ratio between the compact and the porous radii. This quantity is basically equal 
to ^3. Fractal growth is important for small particles creating fluffy agglomerates 
(until 10~^ — 10~^ cm approximately), after this point the relative velocities become 
high enough so compactness becomes important. Before the particles reach a fully 
compacted stage they fragment, become monomers and a new cycle of growth starts. It 
is important to note that the porosity of the aggregates before the first fragmentation 
event is usually higher than the porosity values after equilibrium is reached. This can 
be seen in Figure [2^81 (grains after 400 years and 3000 years). The reason is that before 
the first fragmentation event, particles involved in collisions are typically equal sized 
so these particles produce fluffy structures. However, when the distribution function 
relaxes in equilibrium, there are collisions between smaller and bigger aggregates as well 
which results in somewhat compacted aggregates. 
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2.5.3.2 Model comparison Ormel et al. (120071 ) 



We compare our Monte Carlo code with the one developed by Ormel et al. (|2007l ). 
They use the Minimum Mass Solar Nebula disk model (MMSN) and somewhat different 
dust-parameters which we changed accordingly (distance from the central star = 1 AU, 
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Figure 2.8: This figure shows the radial enlargement of the dust aggregates after 400 
and 3000 years. The x axis is the compact radii of the particles, the y axis is the ratio 
between the compact and the porous radius, this quantity is basicaUy equal to 5*3. 



temperature = 280 K, density of the gas is 8.5 x 10^^*^ g/cm^, gas to dust ratio = 240, a 
= 10~'^; monomer radius = 0.1 /xm, monomer density = 3 g/cm^, surface energy density 
of the monomers = 25 ergs/cm~^). 

They follow particle coagulation at one pressure scale height above the midplane of the 
disk. Because of this if the particles reach a critical stopping time {Train = a/Q, where 
Cl is the Kepler frequency) , the particles rain out meaning that these particles leave the 
volume of the simulation, the distribution function of the dust p articl es is collapsing as 
it can be seen in their figures (Figure 10 and 11 in Ormel et al. (|2007l )). 



We do not include this effect in our model but we stop the simulation at the first rain out 
event and compare our distribution functions until this point. We use 10^ representative 
particles (100 x 100) during the simulation. 

This can be seen at Figure 12.91 Th e read er is advised to examine this figure together 
with Figure 10. c. from Ormel et al. (|2007l ) because this is the figure we reproduced here. 
Furthermore we would like to point out that the scale of the y axis is different in the 
two figures. Our figure shows two orders of magnitude from the normalized distribution 
functions whereas their figure covers more than 10 orders of magnitude from the real 
distribution function. 

Keeping these in mind, one can compare the results of the two Monte Carlo codes. 



The continuous lines at Figure 12.91 in this chapter show the distribution functions at t 
= 10 years (thin line), 100 years (thicker line), 1000 years (thickest line). The dotted 
line shows the distribution function at the time of the first compaction event (t = 1510 
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Figure 2.9: Distribution functions obtained by using Ormel et al. input pa- 

rameters. The continuous hues show the distribution functions at t = 10 years (thin 
hne), 100 years (thicker hue), 1000 years (thickest hue). The dotted hue shows the 
distribution function at the time of the first compaction event (t = 1510 years), the 
dashed hne shows the distribution function at the first rain out event (t = 2900 years). 



years), the dashed Hne shows the distribution function at the first rain out event (t 
2900 years). The same notation is used by Ormel et ah (|2007l ) at Figure 10. c. 



We compared the position of the peaks of the distribution functions and the approximate 
sha pe of the curves. We can conclude that our code reproduces the results of Ormel et 
al. pOOTt ) very well. 



The required CPU time to perform this simulation is only 10 minutes. One might ask 
why the CPU time is almost ten times smaller now? Why do the previous simulations, 
which used the same number of representative particles (10^) and simulated approxi- 
mately the same time interval (3000 years), take so long? The required CPU time does 
not scale linearly with the used number of particles. It scales linearly with the number 
of collisions simulated. The difference between this run an d the previous two simula- 
tions is fragmentation. In the simulations of Ormel et al. (I2OO7I ) no fragmentation is 
happening because the growth timescales are longer. Using our initial parameters, the 
first fragmentation event happens around 1000 years, the number of small particles are 
never completely depleted after this time. As the small particles thereafter are always 
present, the number of collisions will be much higher than before. 

Also note th at the porosities of these particles would be smaller if the model of Ormel 
et al. (j2007l ) included fragmentation (for the reason see Sect. I2.5.3.1|) . 
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Figure 2.10: The evohition of dust particles including the effects of Brownian motion 
and turbulence, porosity and using two different monomer sizes (ai = 0.1/xm and 02 = 
0.4/im). The particle distribution is saved after t = 3x 10° years - dash-dot line, 3 x 10^ 
years - dashed line, 3 x 10^ years - dotted line, and 3 x 10^ years - continuous line. 

2.5.4 Monomer size distribution 

An interesting question which can easily be answered with our method is: How the 
mixture of different sized monomers change the maximum agglomerate size which can 
be reached? As we can see from Equations 12.141 and \2.15\ the rolling energy is lower 
for smaller monomers and of course the number of monomers in an agglomerate is 
much higher if the same agglomerate is built up of lighter monomers. This would mean 
higher fragmentation energy and one would expect that the particles would be harder 
to fragment resulting in bigger grains. 

We performed a simplified simulation to be able to answer this question. Only two 
different monomer sizes are considered here, oi = 0.1/im and 02 = 0.4/im assuming that 
half of the mass (or representative particles) belongs to the small monomers, the other 
half belongs to the big monomers. 

One problem arises here with the rolling energy. The rolling energy changes with 
monomer size, and as our method cannot follow exactly the number of contacts in an 
aggregate and what kind of monomers are connected, we are forced to use an averaged 
rolling energy. One has to carefully consider what the average rolling energy should be. 
In our case, the big monomer is 64 times heavier than the small monomer. Let's assume 
that 50% of the mass of an aggregate is built up from small monomers; on the other 
hand, if we compare the number of different monomers, the small monomers will be 
64 times more numerous than the big ones. This means that the contribution of small 
monomers in the average rolling energy {Erou) should be higher. This can be achieved 
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by using the following weighting: 

Eroll = r Erolli H r Erolh, (2-21) 

ai + 02 ai + 02 

where Eroiii is the rolling energy between monomers with radius oi, EroU2 is the rolling 
energy between monomers with radius 02. 

As we can see in Figure [2.10t the maximum aggregate sizes reached are approximately an 
order of magnitude higher than on Figure 12.71 as it was predicted earlier in this Section. 



2.6 Conclusions and outlook 

We have shown that our representative particle method for aggregation of particles 
in astrophysical settings works well for standard kernels. It has the usual advantages 
of Monte Carlo methods that one can add particle properties easily and without loss 
of computational speed. Moreover, it naturally conserves the number of computational 
elements, so there is no need to "add" or "remove" particles. Each representative particle 
represents a fixed portion of the total mass of solids. 

Our method may have various possible interesting extensions and applications. Here we 
speculate on a few of these. For instance, the fact that each representative particle corre- 
sponds to a fixed amount of solid mass makes the method ideal for implementation into 
spatially resolved models such as hydrodynamic simulations of planetary atmospheres 
or protoplanetary disks. We can then follow the exact motion of each representative 
particle through the possibly turbulent environment, and thereby automatically treat 
the stochastic nature and deviation from a Boltzmann distribution of the motion of 
particles with stopping times of the same order as the turbulent eddy turn-over time. 
It is necessesary, however, to assure that a sufficiently large number of representative 
particles is present in each grid cell of the hydrodynamic simulation. For large scale 
hydrodynamic simulations this may lead to a very large computational demand for the 
coagulation computation, as well as for tracking the exact motion of these particles. If 
strong clumping of the particles happens, however, much of the "action" anyway hap- 
pens in these "clumps" , and it may then not be too critical that other grid cells are not 
sufficiently populated by representative particles. This, however, is something that has 
to be experimented. 

X^ISUSp SSDLU 'LUJOU 

Our representative particle method can in principle also be used to model the sublimation L o 

and condensation of dust grains. If a particle sublimates then the representative particle °- ^ ^ °- 
becomes simply an atom or molecule of the vapor of this process. It will then follow 
the gas motion until the temperature becomes low enough that it can condense again. 
Other representative particles which are still in the solid phase may represent physical 
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particles that can act as a condensation nucleus. Finally, in our method the properties J^^^H*^ e 
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of the particle can not only change due to collisions, but we can easily implement other 
environmental factors in the alteration of particle properties. 

There are two main drawbacks of the method. First, it only works for large particle 
numbers, i.e. it cannot treat problems in which individual particles start dominating 
their immediate environment. Ormel's method and its expected extension do not have 
this problem. Secondly, the method cannot be accelerated using implicit integration, 
while Brauer's method can. 

All in all we believe that this method may have interesting applications in the field of dust 
aggregation and droplet coagulation in protoplanetary disks and planetary atmospheres. 
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Chapter 3 



Mapping the zoo of laboratory 
experiments 



Based on 'The outcome of protoplanetary dust growth: pebbles, boulders, or planetesi- 
mals? I. Mapping the zoo of laboratory collision experiments' by C. Guttler, J. Blum, 
A. Zsom, C. W. Ormel & C. P. Dullemond published in A&A, 513, 56. 



3.1 Introduction 



The first stage of protoplanetary growth is still not fully understood. Although our 
empirical knowledge on the collisional properties of dust aggregates has considerably 
widened over the past years ( Blum Sz Wurm . 20081 ) . there is no self-consistent model for 
the growth of macroscopic dust aggregates in protoplanetary disks (PPDs) . A reason for 
such a lack of understanding is the complexity in the collisional physics of dust aggre- 
gates. Earlier assumptions of perfect sticking have been experimentally proven false for 
most of the size and velocity ranges under consideration. Recent work also showed that 
fragmentation and porosity play important rol es in mutual col li sions between protoplan- 
etary dust aggregates. In their review paper, iBlum fc WurmI (j2008l ) show the complex 
diversity that is inherent to the collisional interaction of dust aggregates consisting of 
micrometer-sized (silicate) particles. This complexity is the reason why the outcome 
of the collisional evolution in PPDs is still unclear and why no 'grand' theory on the 
formation of planetesimals, based on firm physical principles, has so far been developed. 

The theoretical understanding of the physics of dust aggregate collisions has seen major 
progress in recent decades. The behavior of aggregate collisions at low collisional en- 
ergies - where the aggregates show a fractal nature - is theo retically described by the 
molecular dynamics simulations of lDominik Sz Tielena (jl997l ). The predictions of this 
model ~ concerning aggregate sticking, compaction, and catastrophic di sruption - could 
be quantitatively confirmed by the laboratory collision experiments of 
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([2000D- Also, the cohision behavior of macroscopic dust aggregates was successfully 



modele d by a smooth partic le hydrodynamics me t hod , calibrated by laboratory exper- 



iments (iGiittler et alJ . l2009l : iGeretshauser et alJ . 120091 ) . These simulations were able 



to reproduce bouncing collisions, which were observed in many laboratory experiments 



(IBlum fc Wurml . l2008l ). 



As laboratory experiments have shown, collisions between dust aggregates at interme- 
diate energies and sizes are characterized by a plethora of outcomes: ranging from 
(partial) sticking, bouncing, and mass transfer to catastrophic fragmentation fsee Blum 



& Wurm. I2OO8I ). Prom this complexity, it is clear that the construction of a simple the- 
oretical model which agrees with all these observational constraints is very challenging. 
But in order to understand the formation of planetesimals, it is imperative to describe 
the entire phase-space of interest, i.e., to consider a wide range of aggregate masses, 
aggregate porosities, and collision velocities. Likewise, the collisional outcome is a key 
ingredient of any model that computes the time evolution of the dust size distribution. 
These collisional outcomes are mainly determined by the collision velocities of the dust 
aggregates, and these depend on the disk model, i.e. the gas and material density in 
the disk and the degree of turbulence. Thus, the choice of the disk model (including its 
evolution) is another major ingredient for dust evolution models. 

These concerns lay behind the approach we adopt in this and subsequent chapters. That 
is, instead of first 'funneling' the experimental results through a (perhaps ill-conceived) 
theoretical collision model and then to calculate the collisional evolution, we will directly 
use the experimental results as input for the collisional evolution model. The drawback of 
such an approach is of course that experiments on dust aggregate collisions do not cover 
the whole parameter space and therefore need to be extrapolated by orders of magnitude, 
based on simple physical models whose accuracy might be challenged. We still feel that 
this drawback is more than justified by the prospects that our new approach will provide: 
through a direct mapping of the laboratory experiments, collisional evolution models can 
increase enormously in their level of realism. 

In this chapter, we will classify all existing dust- aggregate collision experiments for sili- 
cate dust, including three additional original experiments not published before, according 
to the above parameters (Sect. 13. 2p . We will show that we have to distinguish between 
nine different kinds of collisional outcomes, which we physically describe in Sect. 13. 3i 
For the later use in a growth model, we will sort these into a mass-velocity parameter 
space and find that we have to distinguish between eight regimes of porous and compact 
dust-aggregate projectiles and targets. We will present our collision model in Sect. 13.41 
and the consequences for the porosities of the dust aggregates in Sect. 13.51 In Sect. 13.61 
we conclude our work and give a critical review on our model and the involved necessary 
simplifications and extrapolations. 



In Chapter d] (jZsom et al.l . |2009| ) we will then, based upon the results presented here. 



follow the dust evolution using a recently invented Monte-Carlo approach (Zsom &: 
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Table 3.1: Table of the experiments which are used for the model. 



projectile mass eollision velocity micro- coUisional outcome reference 

nip [g] 11 [cm s-l] gravity (see Fig. I3.lt 



Exp 1 


7.2 ■ 10-^2 _ 7 2 ■ 10-y 


0.1 


- 1 


yes 


SI 


Blum et al. (1998. 20021. 
Wurm & Blum (1998) 


Exp 2 


7.2 ■ 10-12 _ 2.0 ■ 10-10 


10 


- 50 


yes 


SI 


Wurm & Blum (1998) 


Exp 3 


3.5 ■ 10-12 - 3.5 ■ 10-10 


0.02 


- 0.17 


yes 


SI 


Blum et al. (2000). 




1.0 ■ 10-12 - 1.0 ■ 10-10 


0.04 


- 0.46 


yes 


SI 


Krause & Blum (2004) 


Exp 4 


1.2 ■ 10-10 - 4.3 ■ 10-10 


7 - 


1000 


yes 


S2 


Blum & Wurm (2000) 


Exp 5 


2 ■ 10-3 - 7- 10-3 


15 - 


- 390 


yes 


Bl, Fl 


Blum & Miinch (1993) 




10-^- lo-* 


15 - 


- 390 


yes 


Bl, Fl 




Exp 6 


10-" - lo-* 


10 - 


- 170 


yes 


S2, S3 


Lanekowski et al. (2008) 




10-'' - 3 ■ 10-3 


50 - 


- 200 


yes 


B2, S2, S3 






2.5 ■ 10-^ - 3 ■ 10-3 


200 


- 300 


yes 


S3 




Exp 7 


10-3 - 3 . 10-2 


20 - 


- 300 


yes 


S3 


Blum & Wurm (2008) 


Exp 8 


10-3 - 3.2 ■ 10-2 


16 


- 89 


no 


S3 


Giittler et al. (2009) 


Exp 9 


10-3 - 10-2 


10 


- 40 


yes 


Bl 


D. HeiBelmann et al. (in prep.) 




10-3 - 10-2 


5 - 


- 20 


yes 


Bl 




Exp 10 


2 ■ 10-3 - 5 ■ 10-3 


1 - 


- 30 


no 


Bl 


Weidling et aL (2009) 


Exp 11 


1.6 ■ 10-" - 3.4 ■ 10-2 


320 


- 570 


yes 


Fl 


Lammcl (2008J 


Exp 12 


3.5 ■ 10-1^ 


1500 


- 6 000 


no 


F2 


R. Schrapler & J. Blum (in prep.) 


Exp 13 


0.2 - 0.3 


1650 


- 3 750 


no 


F2 


Wurm et al. (2005a) 


Exp 14 


0.2 - 0.3 


350 - 


- 2150 


yes 


F2 


Paraskov et al. (2007) 


Exp 15 


0.39 


600 - 


- 2 400 


no 


S4 


Wurm ct al. (2005b) 


Exp 16 


4 ■ lO-'^ - 5 ■ 10-5 


700 


- 850 


no 


S4 


Teiser & Wurm (2009a) 


Exp 17 


1.6 ■ 10-" - 2.0 ■ 10-2 


100 - 


- 1000 


no 


S4 


Sect. 13.2.2.11 


Exp 18 


10-0 - lO-* 


10 - 


1000 


no 


Bl, S2, S4 


Sect. 13.2.2.21 


Exp 19 


1.5 • 10-3 - 3.2 ■ 10-3 


200 


- 700 


ves 


S4. F3 


Sect. 13.2.2.31 



Dullemond, 120081 ) for three different disk models. This is the first fuhy self-consistent 
growth simulation for PPDs. The results presented in Chapter U] represent the state- 
of-the-art modeling and will give us important insight into questions, such as if the 
meter-size barrier can be overcome and what the maximum dust-aggregate size in PPDs 
is, i.e. whether pebbles, boulders, or planetesimals can be formed. 



3.2 Collision experiments with relevance to planetesimal 
formation 



In the past years, numerous laboratory and space e xperiments on t h e coll isional evo- 
lution of protoplanetary dust have been performed (iBlum &: Wurml . l2008l ). Here, we 
concentrate on the dust evolution around a distance of 1 AU from the solar-type cen- 
tral star where the ambient temperature is such that the dominating material class are 
the silicates. This choice of 1 AU reflects the kind of laboratory experiments that are 
included in this chapter, which were all performed with Si02 grains or other refractory 
materials. The solid material in the outer solar nebula is dominated by ices, which 
possibly have very different material properties than silicates, but only a small fraction 
of laboratory experiments have dealt with these colder (ices, organic materials) or also 
warmer regions (oxides). In Sect. 13.6. 2^ we will discuss the effect that another choice of 
material might potentially have, but as we are far away from even basically comprehend- 
ing the collisional behavior of aggregates consisting of these materials, we concentrate 
in this study on the conditions relevant in the inner solar nebula around 1 AU. 
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SI (hit & stick) 
00 




52 (sticl<ing through surface effects) 

00 

53 (sticl<ing by penetration) 



S4 (mass transfer) 






B1 (bouncing with compaction) 

_o 0-. 

B2 (bouncing with mass transfer) 



F1 (fragmentation) 
O o oO 



F2 (erosion) 



F3 (fragmentation with mass transfer) 

.oO„-- 




Figure 3.1: We classify the variety of laboratory experiments into nine kinds of coUi- 
sional outcomes, involving sticking (S), bouncing (B) and fragmenting (F) collisions. All 
these collisional outcomes have been observed in laboratory experiments, and detailed 
quantities on the outcomes are given in Sect. 



Table [3TT] lists all relevant experiments that address collisions between dust aggregates of 
different masses, mass ratios, and porosities, consisting of micrometer-sized silicate dust 
grains, in the relevant range of collision velocities. Experiments 1-16 are taken from the 
literature (cited in Table EH]), whereas experiments 17 - 19 are new ones not published 
before. In the following two subsections we will first review the previously published 
experiments (Sect. 13.2. ip and then introduce the experimental setup and results of new 
experiments that were performed to explore some regions of interest (Sect. I3.2.2p . All 
these collisions show a diversity of different outcomes for which we classify nine different 
collisional outcomes as displayed in Fig. 13. 1[ Details on these collisional outcomes are 
presented in Sect. 13.31 



3.2.1 A short review on collision experiments 



We briefly review published results of dust-collision experiments here since these deter- 
mine the collisional mapp i ng in Sect. 13.31 and l3.4[ The interested reader is referred to the 
review by lBlum &: WurmI (120081 ) for more information. All experiments are compiled and 
referenced in Table 13.11 where we also list the collision velocities and projectile masses, 
as these will be used in Sect. 13.41 Most of the experiments in Table [XT] (exception: Exp 
10) were performed under low gas pressure conditions to match the situation in PPDs, 
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and most of the experiments were carried out in the absence of gravity (i.e. free falhng 
aggregates or micro-gravity facihties), see Col. 4 of Table [3TT1 For the majority of the 
experiments, spherical monodisperse Si02 monomers with diameters between 1.0 fxm 
and 1.9 fiui were used; some experiments used irregular Si02 grains with a wider size 
distribution centered around ~ 1.0 i-im, and Exp 5 used irregular ZrSi04 with monomer 
diameters in the range 0.2 ... 1.0 /im. 

Exp 1 -4- A well-known growth mechanism for small dust aggregates is the hit-and- 
stick growth, in which the aggregates collide with such a low kinetic energy that they 
stick at each other upon first contact without any restructuring. The first experiments 
to unambiguously show th at the hit-and-sti c k pro c ess is relevant to pr o topla.netary dust 



aggregation wer e those by IWurm BlumI ([199^, iBlum et al.l (Il998l . l200d . [200^) and 
Krause BlumI (120041 ). These proved that, as long as the collision velocities for small 
dust aggregates stay well below 100 cm s^^, sticking collisions lead to the formation of 
fractal ag gregates. This agrees with the molecular-dvnamics simulations bv Dominik 



Wada et al 



(2007 



200J 



sT 2009 ). The various experimental ap- 



k Tielens (|l997l l and 
proaches for Exp 1-3 used all known sources for relative grain velocities in PPDs, 
i.e. Brownian motion (Exp 3), relative sedimentation (Exp 1), and gas turbulence (Exp 
2). In these papers it was also shown that the hit-and-stick growth regime leads to a 
quasi-monodisperse evolution of the mean aggregate masses, depleting small grains effi- 
ciently and rapidly. Fo r collis ions between these fractal aggregates and a solid or dusty 
target, I Blum &: WurmI (j2000l . Exp 4) found growth at even higher velocities, in which 
the aggregates were restru cture. This also agrees with molecular-dynamics simulations 
(jDominik &: Tielend . 119971 ). and so this first stage of protoplanetary dust growth has so 
far been the only one that could be fully modeled. 



Exp 5: iBlum Miinchl (jl993l ) performed collision experiments between free falling 
ZrSi04 aggregates of intermediate porosity (0 = 0.35, where (j) is the volume fraction of 
the solid material) at velocities in the range of 15 - 390 cm s^^. They found no sticking, 
but, depending on the collision velocity, the aggregates bounced {v < 100 cm s^^) or 
fragmented into a power-law size distribution {v > 100 cm s~^). The aggregate masses 
were varied over a wide range (10-^ to 7x10-3 g), and the mass ratio of the two collision 
partners also ranged from 1:1 to 1:66. The major difference to experiments 1-4, which 
inhibited sticking in these collisions, were the aggregate masses and their non- fractal but 
still very porous nature. 

Exp 6-8: A new way of producing highly porous, macr oscopic dust aggregates i cf) = 0.15 
for 1.5 /im diameter Si02 monospheres) as described bv lBlum &: Schrapleii (j2004l ) allowed =0 
new experiments, u sing the 2.5 cm d i amete r aggregates as targets and fragments of ° 
these as projectiles ( Langkowski et al.l. 2008, Exp 6). In their collision experiments in 
the Bremen drop tower, iLangkowski et al.l (|2008l ) found that the projectile may either 
bounce off from the target at intermediate velocities (50 - 250 cm s^^) and aggregate 
sizes (0.5 - 2 mm), or stick to the target for higher or lower velocities and bigger or 
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smaller sizes, respectively. This bouncing went with a previous slight intrusion and a 
mass transfer from the target to the projectile. In the case of small and slow projectiles, 
the projectile stuck to the target, while large and fast projectiles penetrated into the 
target and were geometrically embedded. They also found that the surface roughness 
plays an important role for the sticking efficiency. If a projectile hits into a surface 
depression, it sticks, while it bounces off when hitting a hill with a small radius of 
curvature comparable to that of the projectile. A s i milar behavior for the sticking by 
deep penetration was also found by iBlum &: WurmI (120081 . Exp 7) when the projectile 
aggregate is solid - a mm-sized glass bead in their case. Continuous experiments on 
the penetr ation of a solid projectile ( 1 to 3 mm diameter) into the highly porou s target 
{(p = 0.15, blum fc Schraplerl . boOJ) were performed bv Iciittler et al.l pOod . Exp 8) 
who studied this setup for the calibration of a smoothed particle hydrodynamics (SPH) 
collision model. We will use their measurement of the penetration depth of the projectile. 



Exp 9-10: As a follow-up experiment of the study of iBlum &: MiinchI (|l993l l. D. 
Heifielmann, H.J. Eraser and J. Blum (in prep., Exp 9) used 5 mm cubes of these highly 
porous (0 = 0.15) dust aggregates and crashed them into each other {v = 40 cm s~^) 
or into a compact dust target with <f) = 0.24 {v = 20 cm s~^). In both cases they 
too found bouncing of the aggregates and were able to confirm the low coefficient of 
restitution (fafter/^^before) of £ = 0.2 for central collisions. In their experiments they 
could not see any deformation of the aggregates, due to the limited resolution of their 
camera, which could h ave explained th e dissip ation of energy. This line of experiments 
was taken up again by IWeidling et al.l (|2009l . Exp 10) who studied the compaction of 
the same aggregates which repeatedly collided with a solid target. They found that 
the aggregates decreased in size (without losing significant amounts of mass), which is 
a direct measurement of their porosity. After only 1 000 collisions the aggregates were 
compacted by a factor of two in volume filling factor, and the maximum filling factor 
for the velocity used in their experiments (1 - 30 cm s~^) was found to be = 0.36. In 
four out of 18 experiments, the aggregate broke into several pieces, and they derived a 
fragmentation probability of Pfrag = 10~^ for the aggregate to break in a collision. 



Exp 11: T he same fragments of the high porositv i 
h Schrapler (120041) as well as intermediate porosity 



0.15) dust aggregates of Blum 
= 0.35) aggregat es were used by 



Lammell ()2008l . Exp 11) who continued the fragmentation experiments of 



Blum fc Miinch 



I993I ). For velocities from 320 to 570 cm s ^ he found fragmentation and measured the 
size of the largest fragment as a measure for the fragmentation strength. 

Exp 12 - 14: Exposing the same highly porous (0 = 0.15) dust aggregate to a stream of 
single monomers with a velocity from 1500 to 6 000 cm s~^, R. Schrapler and J. Blum 
(in prep., Exp 12) found a significant erosion of the aggregate. One monomer impact can 
easily kick out tens of monomers for the higher velocities examined. They estimated the 
minimum vel ocity for this process in an analyti cal model to be appro x. 350 cm s^^. On 
a larger scale. IWurm et al.l (j2005al . Exp 13) and lParaskov et al.l (120071 . Exp 14) impacted 
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dust pr ojectiles with masses o f 0.2 to 0.3 g and sohd spheres into loosely packed dust 
targets. 



Paraskov et al 



(|2007l ) were able to measure the mass loss of the target in drop- 
tower experiments which was - velocity dependent - up to 35 projectile masses. The 
lowest velocity in these experiments was 350 cm s~^. 

Exp 15-16: In a collision between a projectile o f intermediate poros ity and a com- 



Wurm et al. 



(|2005bl . Exp 15) found 



pressed dust target at a velocity above 600 cm s~^, 
fragmentation of the projectile but also an accretion of mass onto the target. This ac- 
cretion w as up to 0.6 pro j ectile m asses in a single collision depending on the collision 
velocity. Teiser &: Wurm (j2009al . Exp 16) studied this partial sticking in many colli- 
sions, where solid targets of variable sizes were exposed to 100 to 500 /um diameter dust 
aggregates with a mean velocity of 770 cm s^^. Although they cannot give an accre- 
tion efficiency in a single collision, they found a large amount of mass accretion onto 
the targets, whi ch is a combination of t he pure partial sticking and the effects of the 
Earth's gravity. iTeiser &: WurmI (j2009al ) argue that this acceleration is equivalent to 
the acceleration that micron-sized particles would experience as a result of their erosion 
from a much bigger body which had been (partially) decoupled from the gas motion in 
the solar nebula. 



3.2.2 New experiments 

In this section, we will present new experiments which we performed to explore some 
parameter regions where no published data existed so far. All experiments cover colli- 
sions between porous aggregates with a solid target and were performed with the same 
experimental setup, consisting of a vacuum chamber (less than 0.1 mbar pressure) with 
a dust accelerator for the porous projectiles and an exchangeable target. The accelerator 
comprises a 50 cm long plastic rod with a diameter of 3 cm in a vacuum feed through. 
The pressure difference between the ambient air and the pressure in the vacuum cham- 
ber drives a constant acceleration, leading to a projectile velocity of up to 900 cm s~^, 
at which point the accelerator is abruptly stopped. The porous projectile flies on and 
collides either with a solid glass plate (Sect. 13.2.2.11 and I3.2.2.2p or with a free falling 
glass bead, which is dropped when the projectile is accelerated (Sect. I3.2.2.3p . The 
collision is observed with a high-speed camera to determine aggregate and fragment 
sizes and to distinguish between the collisional outcomes (i.e. sticking, bouncing, and 
fragmentation). The experiments in this section are also listed in Table [3TT] as Exp 17 
to 19. 

3.2.2.1 Fragmentation with mass transfer (Exp 17) 

In this experiment, mm-sized aggregates of different volume filling factors (0 = 0.15 
and = 0.35) collided with a flat and solid glass target and fragmented as the collision 
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Figure 3.2: Example for a collision of a porous (0 — 0.35) aggregate with a solid 
target at a velocity of 620 cm s^^. The aggregate fragments according to a power- law 
size distribution and some mass sticks to the target (bottom frame). 

velocity was above the fragmentation threshold of approx. 100 cm s~^. The projected 
projectile size and its velocity were measured by a high-speed camera (see Fig. 13. 2p . In 
few experiments, the sizes of the produced fragments were measured for those fragments 
that were sharply resolved, which yielded a size distribution of a representative number 
of fragments (the number of resolved fragments varied from 100 to 400). Assuming a 
spherical shape of the fragments and an unchanged porosity from the original projectile, 
we calculated a cumulative mass distribution as shown in Fig. 13. 3( where the cumulative 
mass fraction X^iLo("^i/-^^F) is plotted over the normalized fragment mass my^/nip. Here, 
nii and M-p = X^i^i "^-i mass of the i-th. smallest fragment, and the total mass 

of all visible fragments and is the total number of fragments. We found that the 
cumulative distribution can be well described by a power law 

/ n{m')m! dm' = ( — ] , (3.1) 

where m' and m are the mass of the fragments in units of the projectile mass and ^ 
is a parameter to measure the strength of fragmentation, defined as the mass of the 
largest fragment divided by the mass of the original projectile. The deviation between 
data and power-law for low masses (see Fig. 13. 3p is due to the finite resolution of the 
camera, which could not detect fragments with sizes <C 50 /zm. In the ten experiments 
where the mass distribution was determined, the power-law index k was nearly constant 
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10'^ 10'" 10'^ 10'^ 10'' 10° 
fragment moss / projectile mass 



Figure 3.3: Mass distribution for two experiments at the vefocities of 120 and 640 
cm s^^. For the higher masses, the distribution fohows a power-law, while the lower 
masses are depleted due to the finite camera resolution. The slopes are the same for 
both experiments, and there is only an offset (pre-factor) between the two. The inset 
describes this pre-factor ii (cf. Eq. 13. ip which is a measure for the strength of the 
fragmentation. The value clearly decreases with increasing velocity (Eq. 13. 2|) . 



from 0.64 to 0.93, showing no dependence on the velocity, which varied from 120 to 
840 cm s~^. However, a clear dependence on the velocity was found for the parameter 
/i, which decreased with increasing velocity as shown in the inset of Fig. 13.31 This 
increasing strength of fragmentation can be described as 

a(v) = ( , (3.2) 

' VlOO cm s-i/ ' ^ ^ 

where the exponent has an error of ±0.2. The curve was fitted to agree with the observed 
fragmentation threshold of 100 cm s~^. 

It is important to know that the number density of fragments of a given mass follows 
from Eq. 13.11 as 

n(m') = -^rn'^'^ (3-3) 

and that the power law for this mass distribution can be translated into a power-law size 
distribution n(a) cx with A = 3k — 4. This yields A values from —2.1 to —1.2, which 
is m uch flatter than the power-law index of —3.5 from the MRN distribution ( Mathis 
et ah, Il977l ). which is widely used for the description of high-speed fragmentation of 



solid materials. More over, this power-law index is consistent with measurements of °° 



Blum Sz Miinchl (jl993l ) who studied aggregate-aggregate collisions between millimeter- - 
sized ZrSi04 aggregates (see Sect. 13. 2p . Their power-law index equivalent to A was 
— 1.4, and for different velocities they also found a constant power-law index and a 
velocity-dependent pre-factor (their Fig. 8a). 
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number of impacts 



Figure 3.4: Mass gain of a solid target in 133 collisions (S. Kothe, C. Giittler & J. 
Blum, unpublished data). The target was weighed after every 19 collisions. After 57 
collisions, one projectile mass of dust was chipped off the target, which is a clear effect 
of gravity. Thus, we added this mass to the following measurements (triangles) and 
fitted a linear mass gain, which is 0.023 x nip in every collision (solid line). 



While most of the projectile mass fragmented into a power-law distribution, some mass 
fraction stuck to the target (see bottom frame in Fig. I3.2p . Therefore, the mass of the 
target was weighed before the collision and again after 19 shots on the same spot. The 
mass of each projectile was weighed and yielded a mean value of 3.34 it 0.84 mg per 
projectile. The increasing mass of the target in units of the projectile mass is plotted in 
Fig. 13.41 After 57 collisions, dust chipped off the target, which can clearly be credited 
to the gravitational influence. For the following measurements we therefore added one 
projectile mass to the target because we found good agreement with the previous values 
for this offset. The measurements were linearly fitted and the slope, which determines the 
mass gain in a single collision, was 2.3 % (S. Kothe, C. Guttler & J. Blum, unpublished 
data) . 



3.2.2.2 Impacts of small aggregates (Exp 18) 



Using exactly the same setup as in the previous section, we performed collision experi- 
ments with very small (20 fim to 1.4 mm diameter) but non-fractal projectiles. Those 
aggreg ates were fragments of larger dust samples as described by iBlum &: Schrapler 
(|2004l ) and had a volume filling factor of (p = 0.15. In this experiment we observed not 
only fragmentation but also bouncing and sticking of the projectiles to the solid glass 
target. Thus, the analysis with the high-speed camera involved the measurement of pro- 
jectile size, collision velocity, and collisional outcome, where we distinguished between 
(1) perfect sticking, (2) perfect bouncing without mass transfer, (3) fragmentation with 
partial sticking, and (4) bouncing with partial sticking. The difference between the cases 
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Figure 3.5: Examples for the experimental outcomes in the collisions of small aggre- 
gates with a solid target. The collision can lead to sticking, bouncing, or fragmentation 
(from left to right). The time between two exposures is 2 ms. 




velocity [cm/s] 



Figure 3.6: Overview on collision experiments between 20 to 1400 ^im diameter ag- 
gregates and a solid target, which leads to sticking (diamonds), bouncing (triangles), 
or fragmentation (crosses). The intermediate sticking-bouncing collision is indicated 
by the squared symbols. The color indicates the sticking probability, i.e. the fraction 
of sticking events in a logarithmic bin around every node. The dotted box denotes 
the approximated parameter range and the solid lines denote the threshold between 
sticking, bouncing and fragmentation as also used in Fig. 13.111 



(3) and (4) is that in a fragmentation event at least two rebounding aggregates were 
produced, whereas in the bouncing collision only one aggregate bounced off. 

For the broad parameter range in diameter (20 to 1400 ^m) and velocity (10 to 1 000 cm s^^ 
we performed 403 individual collisions in which we were able to measure size, velocity, 
and collisional outcome. Examples for sticking, bouncing, and fragmentation are shown 
in Fig. 13.51 The full set of data is plotted in Fig. 13.61 where different symbols were used 
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200 300 400 500 600 
velocity [cm/s] 



700 



Figure 3.7: The volume gain of a solid particle colliding with a porous aggregate 
depends on the collision velocity. The data points are mean values of 11, 8, and 7 indi- 
vidual experiments (left to right), thus, the error bars show the la standard deviation 
of velocities and volume gain in these. The images with a width of 1.9 mm show the 
original 1 mm glass bead and examples for the mass gain in the three corresponding 
collision velocities (S. Olliges & J. Blum, unpublished data). 



for different collisional outcomes. Clearly, collisions of large aggregates and high veloc- 
ities lead to fragmentation, while small aggregates tend to bounce off the target. For 
intermediate aggregate mass (i.e. nip = 10"'' g), all kinds of collisions can occur. The 
background color shows a sticking probability, which was calculated as a boxcar average 
(logarithmic box) at every node where an experiment was performed. Blue color denotes 
a poor sticking probability, while a green to yellow color shows a sticking probability of 
approx. 50 %. We draw the solid lines in a polygon [(100,70,800,200,200, 17) cm s~\ 
(1.6 • 10"'', 5 • lO"'^, 1 • 10"'^, 8 • 10~^°, 1 • 10"^ 1 • 10"^) g] to mark the border between 
sticking and non-sticking as we will use it in Sect. 13. 4i For the higher masses, this 
accounts for a bouncing-fragmentation threshold of 100 cm at 1.6 • 10"'* g (Exp 18), 
and for the lower masses, we assume a constant fragmentation threshol d of 200 cm s~^, 
which roughly agrees with the restructuring- fragmentation threshold of iBlum fc Wurm 
([20001, Exp 4). For lower velocities outside the solid-line polygon, bouncing collisions 
are expected, whereas for higher velocities outside the polygon, we expect fragmenta- 
tion. Thus, an island of enhanced sticking probability for 10"^ - 10"'' g aggregates at a 
broad velocity range from 30 to 500 cm s~^ was rather unexpected before. The dotted 
box is just a rough borderline showing the parameters for which the experiments were 
performed as it will also be used in Sect. [37 
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3.2.2.3 Collisions between similar sized solid and porous aggregates (Exp 



In a collision between a free falling glass bead of 1 mm diameter and a porous (0 = 0.15) 
dust aggregate of 1.5 to 8.5 mg mass, we observed fragmentation of the porous aggregate 
while some mass was growing on the solid and indestructible glass bead (S. Olliges & 
J. Blum, unpublished data). In this case, the high-speed camera was used with a 3D 
optics that allowed the imaging of the collision from two angles, separated by 90°. On 
the one hand this made it possible to exactly measure the impact parameter b also if the 
offset of the two collision partners is in the line of sight of one viewing angle. Moreover, 
observing the mass growth of the solid projectile is not only a projection in one direction 
but can be reconstructed to get a 3D measurement. The relative velocity and aggregate 
size were accordingly measured from the images before the collision while the mass gain 
of the solid glass bead was measured after the collision. Figure 13.71 shows a diagram 
of the volume gain in units of projectile volume (projectile: porous aggregate) over the 
collision velocity. The three data points are averaged over a number of experiments at 
the same velocity. The error bars denote the la standard deviation of collision velocities 
and projectile volume, respectively. A clear trend shows that the volume gain of the 
solid particle decreases with velocity, and we fitted the data points with 



where Vp is the volume of the glass bead. In this experiment we were not able to measure 
the size distribution of the fragments because the absolute velocity is determined by the 
projectile velocity (up to 600 cm s~^), and the faster fragments were out of the frame 
before they were clearly separated from each other. 

3.3 Classification of the laboratory experiments 

In this section, the experiments outlined above will be categorized according to their 
physical outcomes in the respective collisions. In Sect. 13.21 we saw that various kinds 
of sticking, bouncing, and fragmentation can occur. Here, we will keep all these experi- 
ments in mind and classify them according to nine kinds of possible collisional outcomes 
that were observed in laboratory experiments. These collisional outcomes are displayed 
in Fig. 13.11 The denomination of the classification follows S for sticking, B for bouncing, 
and F for fragmentation. S and F are meant with respect to the target, i.e. the more % 
massive of the two collision partners. We will discuss each of the pictograms in Fig. 13. H ° 
describe the motivation for the respective collisional outcomes and physically quantify 
the outcome of these collisions. 



19) 




(3.4) 



(1) Sticking collisions: A well-known growth mechanism is d ue to hit-and-stick (S I) col- 
lisions. Hit-and-stick growth was observed in the laboratory (|Blum &: Wurml . 120001 : Blum 
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et al.. l200d ) and numerically described (jPominik &: Tielensl . 119971 ) . Experiments show 
that the mass distribution during the initial growth phase is always quasi-monodisperse. 
The evolution of the mean mass within an ensemble of dust aggregates due to hit-and- 
stick (SI) collision s was calculated t o follow a power-law in t ime, i n good agreement with 



the ex periments (IWurm &: Blurnl. Il998l: iKrause &: Bluml . 12004 ) . 



(|l997l ) showed theoretically and 



Blum &: Wurm ( 



200ol ) confirmed experimentally that 



Dominik &: Tielens 



small fractal aggregates stick at first contact if their collision energy is smaller than a 
threshold energy. For higher energies, experiments s howed that an a, g grega te is elas- 



ticall v and plastically deformed at the contact zone (iBlum Sz Miinchl. Il993l: Weidling 
et al., [20091). This increases the number of contacts, which can then l ead to sticking at 



higher velocities, an effect we call sticking through surface effects (S2). iLangkowski et al 



(|2008l ) found that sticking can occur for even larger velocities if the target aggregate is 
porous and significantly larger than the projectile. In this case, the projectile sticks by 
deep penetration (S3) into the target and cannot rebound simply because of geometrical 
considerations. Thi s effect holds al s o true if the projectile aggreg ate is compact , whic h 



has been shown by lBlum &: WurmI (|2008l ) and further studied by lGiittler et al.l ([20091). 
In Sect. 13.2.2.11 we saw that the growth of a solid target can occur if a porous projectile 
fragments and partia lly sticks to the targ et surface (S4). This grow th mechanism was 
already described by [Wurm et al.[ ([2005bl ). [Teiser &: WurnJ ([2009al ) found it to be an 
efficient growth mechanism in multiple collisions. 

(2) Bouncing collisions: If the collision velocity of two dust aggregates is too low for 
fragmentation and too high for sticking to occur, the dust aggregates will bounce (Bl). 
D. Heifielmann et al. (in prep.) found highly inelastic bouncing between similar-sized 
porous dust aggregates and between a dust aggregate an d a dusty but rather compact 
target, where 95 % of the kinetic energy were dissipated. [Weidling et al.[ ([2009[ ) showed 
that the energy can effectively be dissipated by a significant (and for a single collision 
undetectable) compaction of the porous aggregates after multiple collisions (collisional 
outcome bounc ing with compaction (Bl )). Another kind of bouncing occurred in the 
experiments of [Langkowski et al.[ ([2008[ ) in which a porous projectile collided with a 
significantly bigger and also highly porous target aggregate. If the penetration of the 
aggregate was too shallow for the S3 sticking to occur, the projectile bounced off and 
took away mass from the target aggregate. This bouncing with i nass t ransfer (B2) was 
also observed in the case of compact projectiles ([Blum Sz Wurml . [20081 ). 



(3) Fragmenting collisions: Fragmentation (Fl), i.e. the breakup of the dust aggre- 
gates, occurs in collisions between similar - sized dust aggregates at a velocity above the 
fragmentation threshold. [Blum &: Miinchl ([l993l ) showed that both aggregates are then 
disrupted into a power-law size distribution. If a target aggregate is exposed to impacts 
of single monomer grains or very small dust aggregates, R. Schrapler &: J. Blum (in 
prep.) found that the target aggregate is efficiently eroded (F2) if the impact velocities 
exceed 1 500 cm s~^. This mass loss of the target was also observed in the c ase of larger 



projectiles into porous targets ([Wurm et al 



2005a 



Paraskov et al. 



2003). Similar to 
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the Fl fragmentation, it may occur that one aggregate is porous while the other one 
is compact. In that case, the porous aggregate fragments but cannot destroy the com- 
pact aggregate. The compact aggregate accretes mass from the porous aggregate (Sect. 
I3.2.2.3p . We call this fragmentation with mass transfer (F3). 

These nine fundamental kinds of collisions are all based on firm laboratory results. 
Future experiments will almost certainly modify this picture and potentially add so far 
unknown collisional outcomes to this list. But at the present time this is the complete 
picture of possible collisional outcomes. Below we will quantify the thresholds and 
boundaries between the different collision regimes as well as characterize physically the 
collisional outcomes therein. 



SI: Hit-and-stick growth 



Hit-and-stick growth occurs when the colli s ional energy involved is less than 5 • -EroU 



Dominik &: Tielensl . 



1997 



Blum fc Wurml . l2000l ). where Sroii is the energy which is 



dissipated when one dust grain rolls over another by an angle of 90°. We can calculate 
the upper threshold velocity for the hit-and-stick mechanism of two dust grains by using 
the definition relation between rolling energy and rolling force, i.e. 



TT 



-E-roU — -^O-qF^oH 



(3.5) 



Here, oq is the radius of a dust grain and -FroU is the rolling force. Thus, we are inside 
the hit-and-stick regime if 

(3.6) 



1 2 



where is the reduced mass of the aggregates. The hit-and-stick velocity range is then 
given by 



V < W5 



vrao-Froii 



(3.7) 



S2: Sticking by surface effects 



For velocities exceeding the hit-and-stick threshold velocity (Eq. 13. 7p . we assume sticking 
because of an increased contact area due to surface flattening and, therefore, an increased 
number of sticking grain-grain contacts. For the calcu lation of the contact area, we take 
an elastic deformation of the aggreg ate (|Hertzj . Il88lh and get a radius for the contact 
area of 

' 15' 
32 



So 



G 



(3.8) 



Here, v is the collision velocity, G is the shear modulus, and is the redu ced radius. We 
estimate the shear modulus with the shear strength, which follows after ISirond (j2004l ) 
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as the geometric mean of t he compressive strength and the tensile strength. These pa- 
rameters were measured bv lBlum &: Schraplerl (j2004l ) to be 4000 dyn cm^^ (compressive 
strength) and 10 000 dyn cm~^ (tensile strength), s o we take 6 320 dyn cm~^ for the 
shear modulus, which is consistent with estimates of IWeidling et al.l (|2009l ). 



The energy of a pair of bouncing aggregates after the collision is 

1 



rest. 



2 ^ 2 

£ -rrinV 
2 ^ 



(3.9) 



with the coefficient of restitution e. The contact energy of the flattened surface in contact 
is 

2 

Econt. = So (3.10) 



a, 



2 ' 




where Eq is the sticking energy of a monomer grain with the radius oq- We expect 
sticking for E'cont. > £^rest., thus, 



15 \ m^af^v'^ 



32 



G 



^Eq 2^ 2 

> £ -nif^v or 



a, 







V < 



32 





3 


205^0 


G 







(3.11) 
(3.12) 



This is the sticking threshold velocity for sticking through surface effects (S2), which 
is based on the Hertzian deformation, which is of course a simplified model, but has 
proven as a goo d concept in many attempts to describe sligh t deformation of porous 
dust aggregates (|Langkowski et al" . 2008 : Weidhng et al. . 2009 ). 



We have to ensure that the centrifugal force of two rotating aggregates, sticking like 
above, does not tear them apart, which is the case if 



cent 



(3.13) 



where T is the tensile strength of the aggregate material. The centrifugal force in the 
worst case of a perfectly grazing collision is 



2„,2 



cent 



(3.14) 



where 2a^ is a conservative estimation for the radial distance of the masses with the 
tangential velocity ev. Thus, only collisions with velocities 



V < 



\32 





3 




G 







(3.15) 



can lead to sticking. For the relevant parameter range (see Table [3^ below) . the thresh- 
old velocity in Eq. 13.151 is always significantly greater than the sticking velocity in Eq. 
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\'S.12\ thus, we can take Eq. 13.121 as the relevant velocity for the process S2. 

We will use this kind of sticking not only within the mass and velocity threshold as 
defined by Eq. 13.121 but also for collisions where we see sticking which cannot so far be 
explained by any model, like in experiment 6 or 18. For all these cases, we assume the 
porosity of target and projectile to be unchanged, disregarding any slight compaction as 
needed for the deformation. One exception is the sticking of small, fractal aggregates , 



which clearly goes tog ether with a compaction of the projectile (jPominik &: Tielens , 



1997 



Blum &: Wurml . I2OOO). In these cases we assume a projectile compaction by a factor of 



1.5 in volume filling factor as there is no precise measurement on this compaction. 



S3: Sticking by deep penetration 



If the target aggregate is much larger than the projectile, porous and flat, an impact 
of a (porous or compact) projectile results in its penetration into the target. Sticking 
is inevitable if the penetration of the projectile is deep enough, i.e. deeper than one 
projectile radius. In that case, the projectile canno t bounce off the target fro m geometric 
considerations. This was fo und in experiments of iLangkowski et al.l (j2008l ) in the case 
of porous projectiles and by lBlum &: Wurml ([2003) in the case of solid projectiles. The 
result of the collision for penetration depths Dp > Op is that the mass of the target is 
augmented by the mass of the projectile, and the volume of the new aggregate reads 



V = Vt-TTal{Dp 



ip) + 2^P 



Vt + -Vp - TTa'pDp , 



(3.16) 
(3.17) 



with Vp and Vt being the volume of the projectile and target, respecti vely. We distinguish 
bet ween compact and porous projectiles and take the experiments of 
and ILangkowski et al.l (j2008l ) for impacts into = 0.15 dust aggregates and calculate 



Giittler et al, 



(|2009h 



the sticking threshold velocities. 



For comvact projectiles, we use the linear relation for the penetration depth of Giittler 



et al. (200' 



3 



Dr 



7- 



(3.18) 



4 



where rup = ^7rp ri(j)pap and 
tively. Although iGiittler et al.l (|2009f ) suggest a power-law relation for the penetration 
depth, i.e. Dp = ^77jO-23±o.i3^o.89±o.34^ choose the linear relation in Eq. 13.181 for 
simplicity, which also ag rees with the data w ithin the error bars. For such a linear fit, 
the slope to the data in iGiittler et al.l ( 2009 ) is 7 = 8.3 • 10~^ cm^ s g~^. We assume 



TToi are the projectile mass and cross section, respec- 



sticking for Dp > ap and get sticking due to process S3 in the velocity range 
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which only depends on the projectile bulk density po and filling factor (pp and not on 
the projectile radius. 

A porous projectile, colliding with a porous target, makes a visible indentation into the 
target aggregate if the kinetic energy is E > E'mirn with a material-dependent minimum 
energy -Emin- The crater volume is then given by 

3 

Kr. = (§-] ' cm^ , (3.20) 



(see Fig. 15 in iLangkowski et al.l . l2008l ). Again, from geometrical considerations, we 
assume that sticking occurs if the projectile penetrates at least one radius deep, thus, 
Vcr. > 0.51^, where Vp = fvrap is the volume of the projectile. Thus, 

, 3 

1 

>- E,. (i=)* (3.22) 
-(a)'- 

For these velocities, the projectile is inevitably embedded into the target aggregate. 
However, if the impact energy is less than -Emim the collision will not lead to a penetration 
so that the final condition for sticking of a porous projectile according to process S3 is 



m 'i2p4<^4 



->-axh/^,(:^) 1 . (3.24) 



S4: Partial sticking in fragmentation events 



As introduced in Sect. 13.2.2.11 a fragmenting collision between a porous aggregate and 
a solid target can lead to a partial growth of the target. The mass transfer from the 
projectile to the target is typically 2.3 % of the projectile mass (Fig. 13. 4p . and without 
better knowledge we assume that the transferred mass has a volume filling factor of 
1.5(?!)p. The remaining mass of the projectile fragments according to the power-law mass 
distribution given in Eq. 13.31 with the fragmentation strength from Eq. 13.21 

For a compact projectile aggregate impacting a compact target, the threshold velocity 
for the S4 process is f = 100 cm s~^ and thus identical to that of the Fl process. The 
fragmentation strength is given by Eq. 13.361 
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Bl: Bouncing with compaction 



In a bouncing collision we find compaction of the two collision partners. For similar-sized 



aggregates, the increase of the volume filling factor was formulated by IWeidling et al 
pood , their Eq. 25) to be 



<J)'^{(j),v) > 



(3.25) 



with ^{v) = uq- cm s~^) , (pmaxiv) = 00 + (v/20 cm s"^)"*^^ and uq = 850, 

00 = 0.15, A0 = 0.215 for v < 50 cm s^^. Here, (j)ma.x is the saturation of the filling 
factor aft er many collisions, w hich follows an exponential function with the e-folding 
width (jWeidling et al.l . [20091). In their experiments, v was the velocity of a porous 
projectile colliding with a solid target (infinite mass). In the case of similar-sized colliding 
aggregates, the velocity would be 0.5 • v for each aggregate in a center-of-mass system. 
Therefore, we scale the velocity as 

(3.26) 

nip 

where Vp (vt) is the center-of-mass velocity of the project ile (target). In the case of 
nip <C mt we have the situation of IWeidling et al.l (120091 ) with Vp = v, thus, these 
velocities are chosen to calculate the scaling of u^v) and (/>max(^^) for projectile and 
target compaction, respectively. This means that a projectile with a negligible mass 
with respect to the target cannot compact the target but is only compacted by itself, 
while two aggregates of the same mass are equally compacted. 



For (/!),v.c^(f ). IWeidling et al.l (j2009l ) gave the above relation which is biased by the exper- 
imentally used dust samples and overestimates the compression for very low velocities. 
Therefore, we propose an alternative scaling relation for (pmaxiv)- In a collision with a 
velocity v we can calculate a dynamic pressure 



1 



Pdyn = J^{v) ■ -PV . 



(3.28) 



This pressure is increased bv a factor viv^. as we know from the experiments of Weidling 
et al. ([20091) that the contact area is very small (factor 1/v oi the aggregate surface) 
and that only a very confined volume is compressed. For v = 20 cm s ~^ the pressure 



calculated from Eq. 13.281 is very close to the value given by [Weidling et al.l ([2009[ ). From 
this pressure we cal culate the compression from the compressive strength curve, which 



Giittler et al 



([20091 ) derived for collisions: 



^comp 



>l 



exp 



Igp-lgPn 



+ 1 



(3.29) 
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Figure 3.8: The original corapressive strength curve measured bv lGiittler etaP (|2009l ) 
(Eq. 13.291 sohd hne) is biased by the dust samples used in the experiments. To describe 
also the compression o f dust aggregates with a volume filling factor lower than those 
used bv lGiittler et al. 1 (l2Q09l ). we extrapolate the curve with a power-law (Eq. 13.301 

dashed line) for p < pn^. 



with = 0.12, (j)2 = 0.58, A = 0.58, and pm = 1.3 x 10^ dyn cm~^. This compressive 
strength curve is also biased from the experiments, as its lowest value is (pi = 0.12. 
Assuming the saturation part of the compressive strength curve to be general, we propose 
a power law for p < pm with the same slope as in Eq. 13.29! for (/>comp(Pm)7 which is then 
given by 

'p_ 



i^comp 



4>2 — <l'\ 1 
4>2 + <t>\ ' 2Aln 10 



(3.30) 



and is able to treat the lowest filling factors and pressures. Equations 13.29! and !3.30! 
determine the compression in a confined volume. Taking into account that after many 
collisions only an outer rim of the aggregate is compressed, we reduce the compression 
by a factor /^^ = . 79 to fit the (pma^iv = 20 cm s^^) = 0.365 experimentally measured 



by IWeidling et al. 



Conclusively, we calculate the increase of the volume filling factor from Eq. 13.251 where 
0max is now provided by the dynamical pressure curve as 



</'max(^) — fc ' '/'comp (Pdyn ) 



(3.31) 



where • 



^comp 



is given by Eqs. 13.29! and l3. 301 For the pressure we use Eq. 13.28! and for the 



corresponding velocities we use Eqs. 13.261 and !3.27l to calculate the projectile and target 
compression, respectively. The maximum compression 0niax(^)! which an aggregate can 
achieve in many collisions at a given velocity, is shown in Fig. EI 
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Weidling et al. (|2009l ) found that in this bouncing regime, the aggregates can also frag- 
ment with a low probability. We adopt this fragmentation probability of 



-ffrae 



10" 



(3.32) 



and assume that an aggregate breaks into two similar-sized fragments as suggested by 
their Fig. 5. 



B2: Bouncing with mass transfer 
I ~1 



Langkowski et al. (|2008l ) and lBlum &: Wurml ([200a) found that the collision between a 
projectile (porous or solid) and a porous target aggregate can lead to a slight penetration 



of the projectile into the target followed by the bouncing of th e projectile. This 



a mass transfer from the target to the projectile (see Fig. 7 in iLangkowski et al 



eads to 



20081 



We a ssume that the transferred mass is one projectile mass (Fig. 8 in ILangkowski et al 
20081 ). thus, 



Amt- 



m 



p ' 



(3.33) 



and that the filling factor of the transferred (compacted) material is 1.5 times that of 
the original target material, i.e. 



1.5 X 



(3.34) 



Although the filling factor of the transferred material was not measured, we know that 
the material is si gnificantly compacte d in the collision (see x-ray micro tomography 
(XRT) analysis of lGiittler et al.l . l2009l ). so that the above assumption seems justified. 



Fl: Fragmentation 



When two similar-sized dust aggregates collide at a velocity which is greater than the 
fragmentation velocity of 

100 cm s~\ (3.35) 



rag 



they will both be disrupted. iBlum &: Miinchl (jl993l ) found fragmentation for mm-sized 
ZrSi04 dust aggregates with a porosity of ^ = 0.35 at a velocity greater than 100 cm s~^. 
In their experiments, the aggregates fragmented according to a power-law size distribu- 
tion with an exponent of A = —1.4 (see Sect. 13. 2. 2. ID . which we will use hereafter. The 
two largest fragments togethe r have a mass of iJ ,(v)(m p + nit), where we can determine 
fi{v) from the experirn e nts o f lBlum &: Miinchl (jl993l . ZrSi04 aggregate collisions with 
(j) = 0.35) and iLammell (j2008l . Si02 aggregates of different porosities). These values are 
plotted in Fig. 13.91 and a power-law fit for velocities v > 100 cm s^^ 
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Figure 3.9: The impact strength for aggregate- aggregate cohision also increases for 
higher velocities (decreasing ^, cp. inset in Fig. 13. 3|) . The fitted power- law is given by 

Eq. M 

is shown by the solid line, which is again fitted to match the fragmentation threshold of 
100 cm s~^ (cp. Eq. 13.20 . Here, the error in the exponent is ±0.02. 



F2: Erosion 

If a projectile collides with a significantly larger porous target aggregate at a sufficiently 
high impact velocity, the target may be eroded. R. Schrapler & J. Blum (in prep.) found 
erosion of porous {cp = 0.15) aggregates which were exposed to 1.5 /xm diameter Si02 
monomers (mass mo) at velocities from 1500 to 6000 cm s~^. Their numerical model, 
which fits the experimental data very well, predicts an onset of erosion for a velocity of 
350 cm s^^. The eroded mass grows roughly linear with impact velocity, i.e. 



Am 8 



nip 60 V 100 cm s ^ 



(3.37) 



wher e Am is the amount of eroded mass and rrtp = mn is the projectile mass. Paraskov 



et al. (j2007l ) also found mass loss of a porous target aggregate for velocities from 350 to 
2150 cm s~^, although the process involved is widely different. T hey used porous and 
solid projectiles, and their results (Fig. 4 in iParaskov et al.l . 120071 ) are consistent with 



Am 15 



mp 20 VlOO cm s"! 



(3.38) 



which agrees with non-zero-gravity experiments of IWurm et al.l (l2005al ). who estimated 
a mass loss of 10 projectile masses for velocities of more than 1650 cm s^^. Due to 
the small variation in projectile mass within each of the two experiments, we apply a 
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power-law in mass and merge both experiments to 

(3.39) 



A a , s \ 0.092 

Am D / V \ / nip \ 



nip 80 VlOO cm s ^ ' \mo J 
The velocity range for erosion is therefore 

Ver > 350 cm s"^ (3.40) 

and is consistent in both experiments. 

For compact targets, R. Schrapler & J. Blum (in prep.) were able to measure the velocity 
range for erosion at 

?;er > 2 500 cm s"^ (3.41) 



Due to the nature of the compact target, far less material was eroded, i.e. 

(3.42) 



A c . , / \ 0.092 



rn-p 550 VlOO cm s V^o 

Here, we applied the same power-law index as in Eq. 13.391 due to the absence of large- 
scale experiments in this case. We assume a mass distribution of the eroded material 
according to Eq. 13.21 



F3: Fragmentation with mass transfer 

In Sect. 13.2.2.31 we described the volume transfer from a porous aggregate to a solid 
sphere (assumed to be representative for a compact aggregate) above the fragmentation 
threshold velocity (see Eq. 13. 4p . Without better knowledge, we assume that the trans- 
ferred mass has a volume filling factor of 1.5 times that of the porous collision partner 
(^p) and cannot exceed the mass of the porous aggregate, thus 

Am = m„rt)l-5<?!'p fo.59-6.3 X 10"'' — , (3.43) 
^ ' \ cm s~^ / 

where 'm'p{t) is the mass of the porous aggregate, which can either be projectile or target 

in our definition, depending on its actual mass. For the fragmentation of the porous 

aggregate we assume a power-law distribution following the Fl case. If the collision 

velocity is higher than 940 cm s~^, Eq. 13.431 yields no mass gain for the compact 

aggregate, thus, the mass of the compact aggregate is conserved and only the porous xtisusp ssdlu mjou 

aggregate fragments. ° ° ° ° ° 

p IT) q 

d d d 
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Figure 3.10: Our model distinguishes between porous and compact aggregates, which 
leads to the displayed four types of collisions ( 'pp ', 'pc ', 'cp ', 'cc ') if the collision partners 
are not too different in size (left). The size ratio of projectile and target aggregate was 
identified as another important parameter and we distinguish between similar-sized 
and different-sized collision partners. Thus, in addition to the four collision types on 
the left, impacts of projectiles into much larger targets {'pP', 'pC\ 'cP\ 'cC; the 
target characterized by a capital letter) can also occur (right). The boundary between 
similar-sized and different-sized aggregates is given by the critical mass-ratio parameter 
Tm. Collisions on the left are restricted to nip < mt ^ ''m'Tip, collisions on the right 

happen for mt > Vn^mp. 



3.4 Collision regimes 



In this section we intend to build on the physical descriptions, which we have derived 
in the previous section, and develop a complete collision model for the determination 
of the collisional outcome in protoplanetary dust interactions (Fig. 13. ip . This means 
that for each collision which may occur, a set of collision parameters will be provided 
as input for a numerical model of the evolution of protoplanetary dust (see Chapter 
U]). The most crucial parameters that mainly determine the fate of the colliding dust 
aggregates in each collision are the respective dust-aggregate masses and their relative 
velocity. 

Moreover, in Sect. 13.21 and 13.31 we saw that the porosity difference between the two 
collision partners also has a big impact on the collisional outcome. The only difference 
between the outcomes Fl and F3 (and between S3 and S4) is that the target aggregate is 
either porous or compact. Thus, we define a critical porosity 0c to distinguish between 
porous or compact aggregates. This valu e can only roughly be con fined between (p = 0.15 
(S3 sticking, clearly an effect o f porosity. iLangkowski et al.l . 120081 ) and (p = 0.64 (random 
close packing, clearly compact iTorquato et al.l . l200d ) , and without better knowledge we 
will choose (he = 0.4. 



Another important parameter is the mass ratio of the collision partners. Again, the 
sticking by deep penetration (S3) occurs for the same set of parameters as the fragmen- 
tation (Fl), and only t he cr itical mass ratio Vm = m t/rrip is different. From the work of 
Blum &: Miinchl (|l993l ) and ILangkowski et al.l (|2008l ). we can confine this parameter to 
the range 10 < rm < 1 000 and will also treat it in Chapter [H as a free parameter (with 
fixed values = 10, 100, 1 000). 
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A further parameter, which has an impact on the coUisional outcome, is the impact angle, 
but at this stage we will treat all collisions as central collisions due to a lack of information 
of the actual influence o f the impact angle on th e col lisional r e sult. Experiments by 



Blum Sz Miinchl (119931 ). iLangkowski et al.l (120081 ) . or iLammell (j2008l ) indicate rather 



small differences between central and grazing collisions, so that we feel confident that 
the error due to this simplification is small. Another parameter, which we also neglect at 



this point due to a 



ack o f experimental data, is the surface roughness of the aggregates. 



Langkowski et al.l (120081 ) showed its relative importance, but a quantitative treatment 
of the surface roughness is currently not possible. 

The binary treatment of the parameters (pc and rm leads to Fig. 13.101 whereafter we 
have four different porous-compact combinations and, if we take into account that the 
collision partners can either be similar-sized or different-sized, we have a total of eight 
collision combinations. We will call these 'pp', 'pP', 'cc', 'cC, 'cp', 'cP', 'pc', and 'pC. 
Here, the first small letter denotes the porosity of the projectile {'p' for porous and 
'c' for compact) and the second letter denotes the target porosity, which can be either 
similar-sized (small letter) or different-sized (capital letter). Aggregates with porosities 
(j) < (pc are 'porous ', those with (p > (pc are 'compact '. If the mass of the target aggregate 
iTT-t < ^m^Wp, we treat the collisions as equal-sized, for mt > Tmmp, the collisions are 
treated as different-sized. 

For each combination depicted in Fig. 13.101 we have the most important parameters 
(1) projectile mass rup and (2) collision velocity v, which then determine the collisional 
outcome. As shown in Fig. 13.111 we treat each combination from Fig. 13.101 separately 
and define the collisional outcome as a function of projectile mass and collision velocity. 
For the threshold lines and the quantitative collisional outcomes we use a set of equations, 
which were given in Sect. 13.31 For a quantitative analysis and application to PPDs (see 
Chapter d]) , knowledge of the material parameters of the monomer dust grains and dust 
aggregates is required. In Table 13.21 we list all relevant parameters for 1.5 fim Si02 
spheres, for which most experimental data are available. However, we believe that the 
data in Table 13.21 are also relevant for most types of micrometer-sized silicate particles. 

The only collisional outcome, which is the same in all regimes, is the hit-and-stick (SI) 
process, which, due to its nature, does not depend on porosity or mass ratio but only on 
mass and collision velocity. Thus, all collision combinations in Fig. 13.111 have the same 
region of sticking behavior for a mass-velocity combination smaller than defined by Eq. 
13.71 This parameter region is marked in green because hit-and-stick (SI) can in principle 
lead to the formation of arbitrary large aggregates. Marked in yellow are collisional °o 
outcomes, which do not lead to further growth of the target aggregate, but conserve the ° 
mass of the target aggregate, which is only the case for bouncing with compaction (Bl). 
For simplicity, the weak fragmentation probability of -Pfrag = 10~^ (see Sect. 13. 3|) has 
been neglected in the coloring. The red-marked regions are parameter sets for which the 
target aggregate loses mass. 
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10"* 10"^ 1C° 10^ 10"" 10"^ •. O" 10^ 

velocity [cm/s] velocity [cm/s] 



Figure 3.11: The resulting collision model as described in this chapter. We distinguish 
between similar-sized (left column) and different-sized (right column) collision partners, 
which are either porous or compact (also see Fig. I3.10p . For each case, the important 
parameters to determine the coUisional outcome are the projectile mass and the collision 
velocity. Collisions within green regions can lead to the formation of larger bodies, while 
red regions denote mass loss. Yellow regions are neutral in terms of growth. The dashed 
and dotted boxes show where experiments directly support this model. 

The dashed and dotted boxes m Fig. 13.111 mark the mass and velocity ranges of the 
experiments from Table IXTl In Chapter HJ this plot will help us to see in which parameter 
regions collisions occur and how well they are supported by experiments. We will now 
go through all of the eight plots in Fig. 13.111 and explain the choice for the thresholds 
between the collisional outcomes. 

'pp': In addition to the omnipresent hit-and-stick (SI) regime, which is backed by 
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Table 3.2: Particle and aggregate material properties used for generating Fig. 13.111 



symbol value 



reference 



monomer-grain properties: 

ao 0.75 /um 

mo 3.18 X 10-12 g 

po 2 g cm^^ 

Eo 2.2 X IQ-s erg 

Froii 10-^ dyn 



Blum k Wurm tOOd ). 



iPoppe et al.1 (1200' 




Heim et al.' (|l999h 



aggregate properties: 
? 0l)5 



T 



7 
A 

/c 
1^0 



6320 dyn cm-^ 
10^ dyn cm~2 



0.40 

10 - 1 000 
8.3 X 10^3 



s cm g 



3.5 X 10^ erg 

3.1 X 10-2 erg 

0.12 
0.58 
0.58 

1.3 X 10^ dyn cm 

0.79 

850 

-1.4 



Blum S Munch 
(,199,^ . D. HeiBelmann 
et al. (in prep.) 
this work 



Blum S Schrapler 
(|2004h 
this work 
this work 

'uuttler et aE* (j2009h 

' Langk owski et al.' 
(|2008h 

Langkowski et al.' 
(|2008h 

'uuttler et a? (|2009h 
'uuttler et a? (j2009h 
'uuttler et a? (|2009h 

'uuttler et a? (|2009h 

this work 

'Weidlmg et 5? toO^h 
this work 



experiments 1 ~ 3 in Table ISTTl collisions of porous projectiles can also lead to sticking 
through surface effects (S2), whose threshold is determined by Eq. I3.12[ For higher 
velocities {v > 100 cm s^^, Eq. I3.35p . fragmentation sets in. Bouncing (Bl) and 
fragmentation (Fl) in this regime are well-tested by experiments 5, 9, and 11 in Table 



'pP': As the projectiles are also porous here, we have the same sticking through surface 
effects (S2) threshold as in 'pp'. The same collisional outcome (but with compaction 
of the projectile) was found for collisions of srn a. 



experiment 4 in Table 13. ip . iLangkowski et al.l (|2008l ) (experiment 6) found the S2 



1 agg regates (jBlum Wurml . I200C , 
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collisional outcome for projectile masses 

mp < 10"^ g , (3.44) 

thus we have a horizontal upper limit for S2 in the 'pp ' plot of Fig. 13.111 Extrapolation 
of Exp. 6 to large aggregate masses 

mp > 0.1 g (3.45) 

results in bouncing with mass transfer (B2). A linear interpolation between perfect 
sticking for rup < 10"^ g and p erfect bouncing for rwp > 0.1 g, justified by the sticking 



probabilities shown in Fig. 5 of iLangkowski et al.l (|2008l ). gives a sticking probability for 



the mass range 10 ^ g < rUp < 0.1 g (striped region in the 'pP' of Fig. 13. lip of 

P.tick = -:jlogio TTT^ . (3.46) 



3 VO-1 g 

In Sect. 13.31 we defined the threshold for sticking by deep penetration (S3) by Eqs. 
13.231 and 13. 24t which are prominent in the 'pP' plot for high velocities. For even higher 
velocities, we have erosion of the porous aggregate (F2), defined by the threshold velocity 
in Eq. 13.391 and based on experiments 12 - 14 in Table 13.11 

'cc': Our knowled ge about collisions betw een similar-sized, compact dust aggregates 
is rather limited. Blum &: Miinch (j 19931 ) performed collisions between similar-sized 
aggregates with (f) = 0.35. Although this is lower than the critical volume filling factor 
<j)c as defined in Table [3^21 we assume a similar behavior also for aggregates with higher 
porosity. Therefore, without better knowledge, we define a fragmentation threshold as 
in the 'pp' regime, and take the hit-and-stick (SI) threshold for low energies. We omit 
the sticking through surface effects (82) in this regime because of the significantly lower 
compressibility of the compact aggregates. 

'cC: In this collision regime the experimental background is also very limited. For low 
collision energies we assume a hit-and-stick (SI) growth, for higher velocities bouncing 
with compaction (Bl) and, if the fragmentation threshold [v > 100 cm s^^, Eq. 13.350 is 
exceeded, fragmentation with mass transfer (S4). Based on experiment 12, we have an 
erosion (F2) limit for velocities higher than 2 500 cm s^^ (Eq. I3.4ip . 

'cp'and 'pc': These two cases are almost identical, with the only difference that the 
compact aggregate can either be the projectile or the target (i.e. slightly lower or higher 
in mass than the target aggregate) . The mass ratio of both aggregates is however within 
the critical mass ratio rm. Besides the already-discussed cases SI, 82, and Bl, we assume 
fragmentation 100 cm (Eq. I3.35p . Due to the nature of the collision between a 
compact and a porous aggregate, only the porous aggregate is able to fragment, whereas 
the compact aggregate stays intact. If the compact aggregate is the projectile, the target 
mass is always reduced, thus we have fragmentation with mass transfer (F3) from the 
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target to the projectile. If the target is compact, it grows by fragmentation with mass 
transfer (S4) if the velocity is less than 940 cm s~^ (see Eq. I3.43p . For higher velocities, 
Eq. 13.431 yields no mass gain and so this region is neutral in terms of growth. Collisions 
at high velocities are confirmed by Exp. 19 in this regime. 

'cP': While small collision energies lead t o hit-and-stick (SI), higher energies result in 
bouncing with mass transfer (B2) (Exp. S. lBlum &: Wurml . l2008l ). This region is confined 



by the sticki ng by deep penetrati on (S3) threshold velocity as defined in Eq. 13.191 based 
on Exp. 7 ( Giittler et al.l . l2009l ). At even higher velocities of above 350 cm s~^ (Eq. 
I3.39p . we get erosion of the target aggregate as seen in Exp. 12 - 14. 



'pC: This plot in Fig. 13.111 looks the most complicated but it is supported by a large 
number of experiments. For low collision velocities, we again have hit-and-stick (SI) and 
sticking through surface effects (S2)as well as a transition to bouncing with compaction 
(Bl) for larger collision energies. The existence of the B l bouncing regi o n has been 
shown in Exp. 9 and 10 (D. HeiBelmann et al., in prep.; IWeidling et al.l . [20091). For 
higher velocities and masses above 1.6 • 10"'' g we assume a fragmentation threshold of 
100 cm s^^ with a mass transfer to the target (S4), as seen in Exp. 16 (Sect. I3.2.2.ip . 
For lower masses, the odd-shaped box of Exp. 18 is a direct input from Sect. 13.2.2.21 (see 
Fig. 13. 6p . In the striped region between Bl and S4, we found a sticking probability in 
Exp. 18 of Pstick = 0.5. For lower masses, Exp. 4 showed sticking through surface effects 
(S2) with a restructuring (compaction) of the projectile. As in the 'pP' regime, we set 
the threshold for a maximum mass to 8 • 10"^° g, while t he upper velocity thre shold 
- which must be a transition to a fragmentation regime (|Blum Wurml . |2000| ) - is 
200 cm s~^ from Exp. 4 and 18. 



3.5 Porosity evolution of the aggregates 



Since the porosity of dust aggr egates is a key factor for the outcome of dust aggregate 
collisions (IBlum &: Wurml. I2008D . it is paramount that collisional evolution models follow 



its evolution (Ormeletal 



20071 . Paper II). Therefore we want to concentrate on the 



evolution of the dust aggregates' porosities and recapitulate the porosity recipe as used 
in Sect. 13.31 In this chapter we have used the volume filling factor </) as a quantitative 
value, defined as the volume fraction of material (one minus porosity). In Chapt er [H we 
will also use the enlargement parameter ^ as introduced by 
is the reciprocal quantity = 



Ormeletal 



(l2007l ). which 



Starting the growth with solid dust grains, we have a volume filling factor of 1, which - 
will however rapidly decrease due to the hit-and-stick (SI) growth, prod u cing highly 
porous, fractal aggregates. Here, we use the porosity recipe of Ormel et al. (j2007l ). who 
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describe this fractal growth by their enlargement parameter as 

= ^^,%±^ X f 1 + ^) + . (347) 



where ^'add is a correction factor in case of nip ~ mo and otherwise zero (for details 
see their Sect. 2.4). This equation predicts an increasing porosity in every hit-and-stick 
(SI) collision. In collisions that lead to sticking through surface effects (S2), we assume 
that the compaction of the aggregates is so small, that their porosity is unaffected. So 
the aggregates are merged and only the mass and volume of both are being added, thus, 

^„ = Wi±Ml . (348) 

n + 

One exception for the sticking through surface effects (S2) occurs in a small parameter 



space which is determined by the experim ents oflBlum &: WurrnI (|200d ) . For the smallest 



masses and a velocity around 100 cm s ^ . iBlum &: WurmI (2000) found sticking of fractal 



aggregates in the 'pP' and 'pC regimes, which goes with a restructuring and, thus, 
compaction of the projectiles. In this case, we assume a compaction of the projectile by 
a factor of 1.5 in volume filling factor, thus 

, _ Vt4>t + min (l.SFp^p, (t>c) , . 

An increasing filling factor is also applied for sticking by deep penetration (S3). Here, 
the mass of the projectile is added to the target while the new volume must be less than 
T4 + The new volume filling factor will be 



Vt<j)t + Vp(l)p 



(3.50) 



where Vnew is taken from Eq. 13.171 (compact projectile) or as Vncw = Vt — Vcr. with 
Vcr. from Eq. 13.201 (porous projectile). In the cases where we transfer mass from one 
aggregate to the other, we always assume that this mass is previously compacted by a 
factor of 1.5 in volume filling factor, but cannot be compacted to more than the critical 
filling factor (pQ. For the bouncing with mass transfer (B2) we have good arguments 
for t his assumption as this compaction is consistent with XRT measurements of Guttler 



et al. (|2009[ ). who also showed that it is likely that this compacted material is transferred 



to the projectile (see their Figs. 7 and 9). Without better knowledge, we assume the 
same compaction of transferred material for fragmentation with mass transfer (F3 and 
S4) , and for these three cases we again use Eq. 13.491 It is necessary to swap the indices of 
target and projectile in the case of bouncing with mass transfer (B2) and fragmentation 
with mass transfer (F3), as the projectile is accreting mass in this collisional outcome. 
For the fragments in S4 and F3 as well as for those in the case of Fl and F2, we assume 
an unchanged porosity with respect to the destroyed aggregate. The most sophisticated 
compaction model is used for collisions that lead to bouncing with compaction (Bl). 
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Table 3.3: Overview of the porosity evolution in the different collisional outcomes. 



collisional outcomes 


porosity evolution 


equation 


Q1 
oi 


nuniei 


io.4n 


Co 
DZ 


neutral or compaction 


13.481 or 13.491 




compaction 


13.171 13.2UI I3.5(JI 


o4 (^target ) 


fluffier 


13.491 


( TTXTCM (Ik /■ ' 1 1 

0'4: I piUJCCLilt: ) 


llcULiai 




Bl 


compaction 


13.251 -13.311 


B2 (target) 


neutral 




B2 (projectile) 


both 


|3.49f 


Fl 


neutral 




F2 


neutral 




F3 (target) 


fluffier 


|3.49f 


F3 (projectile) 


neutral 





^The indices of target and projectile must be swapped here. 



Although IWeidling et al.l (120091 ) measured the compaction only for a small range of 
aggregate sizes and collision velocities, they derived an analytic model to scale this 
compaction in collision velocity and showed that it is independent in aggregate mass. 
We follow this model but release it from the experimental bias due to the (j) = 0.15 
samples they used. As outlined in detail in Sect. 13.31 we basically use Eq. 13.251 and 
scale the (/'max(^') according to Eq. 13.311 (furthermore using Eqs. 13.261 - [3T30|) . 

In summary, one can say that the aggregates' porosities can only be increased by the 
collisional outcomes SI, S4, and F3 (see Table 1331) . where the hit-and-stick (SI) collisions 
will have the most effect. While some collisional outcomes are neutral in terms of porosity 
evolution (Fl and F2), the main processes which lead to more compact aggregates are 
S3 and Bl. 



3.6 Discussion 



In the previous sections we have developed a comprehensive model for the collisional 
interaction between protoplanetary dust aggregates. The culmination of this effort is 
Fig. 13.111 which presents a general collision model based on 19 different dust-collision 
experiments, which will be the basis for Chapter [H Since it plays a vital role, it is worth 
a critical appraisal. We want to discuss the main simplifications and shortcomings of 
our current model in a few examples. 

^ A^ISUSp SSDLU 'LUJOU 

"00000 

(1) The categorization into collisions between similar-sized and different-sized dust ag- 5 A S ^ S 



gregates (see Figs. 13.101 and l3.1ip is well-motivated as we pointed out in Sect. 13.41 Still 
we may ask ourselves whether this binarization is fundamentally correct if we need more 
than two categories, or 'soft' transitions between the regimes. At this stage, a more 
complex treatment would be impractical due to the lack of experiments treating this 
problem. 
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(2) The binary treatment of porosity (i.e. (p < (pc for 'porous' and (p > (pc for 'compact' 
dust aggregates) is also a questionable assumption. Although we see fundamental dif- 
ferences in the collision behavior when we use either porous or compact targets, there 
might be a smooth transition from the more 'porous' to the more 'compact' collisions. In 
addition to that, the assumed value (pc = 0.4 is reasonable but not empirically affirmed. 
On top of that, the maximum compaction that a dust aggregate can achieve in a collision 
depends on many p arameters, such as, e.g., the size distribution of the monomer grains 



Blum et al 



20061) and the ability of the granular material to creep sideways inside a 



dust aggregate (jGiittler et al.l . l2009f ) . 



(3) Although the total number of experiments upon which our model is based is unsur- 
passedly large, the total coverage of parameter space is still small (see the experiment 
boxes in Fig. 13. lip . Thus, we sometimes apply extrapolations into extremely remote 
parameter-space regions. Although not quantifiable, it must be clear that the error of 
each extrapolation grows with the distance to the experimentally confirmed domains (i.e. 
the boxes in Fig. 13. lip . Clearly, more experiments are required to fill the parameter 
space, and the identification of the key regions in the mass-velocity plane is exactly one 
of the goals of Chapter [H 

(4) With such new experiments, performed at the 'hot spots' predicted in Chapter [H 
we will not only close gaps in our knowledge of the collision physics of dust aggregates 
but will most certainly reveal completely new effects. That the 'cc' panel in Fig. 13.111 
is rather simple compared to the more complex 'pC is due to the fact that there are 
hardly any experiments that back-up the 'cc' regime, whereas in the 'pC case we have 
a pretty good experimental coverage of the parameter space. 

In summary, the sophisticated nature of our collision model is both its strength and 
its weakness. The drawbacks of identifying four parameters that shape the collision 
outcome are that rather crude approximations and extrapolations have to be made. But 
it is still preferable to acknowledge the role of, e.g., porosity through a binary treatment 
than not to treat this parameter at all. Our new collision model represents the first 
attempt to include all existing laboratory experiments (for the material properties of 
interest); collisional evolution models can enormously profit from this effort. 



3.6.1 The bottleneck for protoplanetary dust growth 

We have presented the framework and physical background for an extended growth 
simulation. What is to be expected from this? This is the place to speculate under 
which conditions growth in PPDs is most favorable. A look at Fig. 13.111 immediately 
shows that large dust aggregates can preferentially grow for realistic collision velocities 
in the 'cC and 'pC collision regimes (and to a lesser extent in the 'pc' case), due 
to fragmentation with mass transfer (S4). A broad mass distribution of protoplanetary 
dust must be present to make this possible. This prerequisite for efficient growth towards 
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planetesimal sizes has also been suggested bv lTeiser &: WurmI (j2009bl . see their Fig. 11). 
Agglomeration experiments with micrometer-sized dust grains and a sticking probability 
of unity (experiments 1 - 3 in Table [3TT]) have shown that nature chooses a rather narrow 
size distribution for the initial fractal growth phase. To see if this changes when the 
physical conditions leave no room for growth under quasi-monodisperse conditions, i.e. 
whether nature is so 'adaptive' and 'target-oriented' to find out that growth can only 
proceed with a wide size distribution, will be the subject of Chapter IU in which we 
apply the findings of this chapter to a collisional evolution model. 



3.6.2 Influence of the adopted material properties 



The choice of material in our model is 1.5 i-im diameter silica dust, as most of the un- 
derlvi ng experiments were performed with this material. Manv experiments f Blum 
Wurm, 



2000 



Langkowski et al. 



2008 



Blum &: Wurm 



20081) showed that this material is 



at least in a qualitative sense representative for other silicatic materials - also for irregu- 
lar grains with a broader size distribution. Still, the grain size of the dust material may 
have a quantitative influence on the collisional outcomes. For example, dust aggregates 
consi sting of 0.1 /xm are assumed to be stickier and more rigid (IWada et al.l . 120071 . 12008 . 
20091 ) ■ because the grain size may scale the rolling force or breaking energy entering 
into Eqs. 13.71 and I3.12[ However, due to a lack of experiments with smaller monomer 
sizes, we cannot give a scaling for our model for smaller monomer sizes at this point. 
Moreover, organic or icy material in the outer regions of PPDs or oxides and sintered 
material in the inner regions may have a big impact on the collisional outcome, i.e. in 
enhancing the stickiness of the material and thereby potentially opening new growth 
channels. 



As for organic materials, 



Kouchi et al 



([2002) found an enhanced sticking of cm-sized 
bodies covered with a 1 mm thick layer of organic material at velocities as high as 
500 cm s~^ and at a temperature of ~ 250 K. Icy materials are also believed to have 

U2m 



Hatzes et al 



an enhanced sticking efficiency compared to silicatic materials 
collided 5 cm diameter solid ice spheres, which were covered with a 10 - 100 //m thick 
layer of frost. They found sticking for a velocity of 0.03 cm s~^, which is in a regime 
where our model for refractory silicatic material predicts bouncing (see 'pp' or 'cc' in 
Fig. 13. lip . Sintering of porous dust aggregate may occur in the inner r egions near the 
centr al star or - triggered by transient heating events (e.g. lightning, 



Guttler et al. 



3,|20q3) 



20081 ) - even further out. Ongoing studies with sintered dust aggregates (jPoppd . 
show an increased material strength (e.g. tensile strength) by an order of magnitude (C. 
Giittler & J. Blum, unpublished data). This would at least make the material robust 
against fragmentation processes and qualitatively shift them from the porous to the 
compact regime in our model - without necessarily being compact. Due to a severe lack 
of experimental data for all these materials, it is necessary and justified to restrict our 
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model to silicates at around 1 AU, while it is to be kept in mind that these examples of 
rather unknown materials might potentially favor growth in other regions in PPDs. 
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Chapter 4 



Introducing the bouncing barrier 



Based on 'The outcome of protoplanetary dust growth: pebbles, boulders, or planetesi- 
mals? II. Introducing the bouncing barrier' by A. Zsom, C. W. Ormel, C. Guttler, J. 
Blum k C. P. Dullemond published in A&A, 513, 57. 



4.1 Introduction 



(Mizuno. 


1980; 


Pollack et al.. 


1996) 



planets are the outcome of an accretion process that starts with micron-size dust grains 
and covers 40 magnitudes in mass. The paradigm can be divided into three stages. 
The first stage involves the formation of rocky planets and the rocky cores of gas giant 
planets and begins with the co agulation of dust in the protoplanetary disks sur rounding 



many pre-m ain-seauence stars (ISafronov 



Il96fl 



Weidenschilling fc: Cuzzi 



1993: Blum k 



Wurm, 120081 ). The next stage of planet formation is the formation of protoplanetary 
cores from the planetesimals. The idea is that the kilometer-size planetesimals are so 
large that gravity begins to dominate leading to the gr avitationa l agglo meration of these 
bodies to rocky planets. This scenario was studi e d by 



U9m 



Schmitt et al 



Lin (|2004h 



i ng nu r nerical method s by Weidenschilling J 198C ) . iNakagawa et al 



Tanaka et al 



Safronov ( 1969 ) and modeled us- 



(119971). IWetherilll (Il990h. iNomura k Nakagawal (120061). Garaud k 



(1983 ).lMizuno et al 



(|2005l ). and several additional authors. These models solve for 
the size distribution of dust aggregates in the disk as a function of time, and investi- 
gate if, where, and how larger dusty bodies form, and how long this takes. Finally, in 
the third stage, gas accretes onto these protoplanets forming giant planets or - in the 



chaotic, gian 


t imoact ohase. until orbital stabilitv is achieved ( 


Chambers . 


2001: 


et al.. 


2006|; 


Thommes et al.. 


2008). 



Kokubo 
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In this study, we focus on the first phase and study the effectiveness of the dust growth 
by surface force, that is, how large particles become due to simple sticking processes 
only. It is known that initially, for micron-size grains, the growth is driven by Brownian 
motio n. This typically leads to slow c ollisio ns and forms aggregates of fractal struc- 



ture (IKempf et al. 



1998 



Blum et al 



19961 ). In the current picture of dust growth. 



as these aggregates grow, at some point the growth will leave the fractal regime, and 
collisi oiis will start to lead to the compaction and breaking of the aggregates f Blum 



Wurm, bood)^ embedding of small bodies into larger aggregates (leading to 'filling up' 
of th ese larger aggregates and compaction caused bv the force of the collision f Ormel 



( Safronov 



et al.. 120071). A s the size of the dust a-ggreg ates increases, di fferential vertical settling 

and turbulence ( Volk et al. 5^ 198ol: Mizuno 



). A s tne size or tne dust aggreg ; 
19691). radial drift JwhiDplel.ll972l) 



et al.. Il988l : lormel fc Cuzzi I2OO7I ) will become important new mechanisms driving rel- 
ative velocities between aggregates. The increasing relative velocities caused by these 
mechanisms will at least partly compensate the lower collision probability due to lower 
surface-over- mass ratio of large aggregates. When the aggregates grow to sizes of be- 
tween mi llimeter and meter, however, the sticking efficiencv drops stronglv (e.s.. Blum 



& Mii nch (119931)) and t he re lative velocities become so large that aggregates can frag- 
ment (jBlum &: WurmI (|2008l ). so-called 'fragmentatio n barrier'). A nother hurdle that 



the particles have to circumvent is the 'drift barrier' (jWeidenschillingl . Il977al ). namely 
that millimeter or centimeter-siz ed particles are lost to the star because of radial drift 
that occurs on a short timescale. lOkuzumil (j2009l ) pointed out the existence of a 'charge 



barrier', which possibly halts the particle growth at an early stage of fractal aggregates. 
Despite many years of efforts, it is unknown whether the coagulation process can over- 
come these barriers. These barriers have been and remain the main open question about 
the initial stages of planet formation, i.e., the growth from dust to planetesimals. 



Several mechanisms have been pr oposed to overcome this problem, among wh i ch are 



the trapping of dust in vortices (jBarge &: Sommerial . Il995l : iKlahr &: Henning . 



1997 



Lvra et al.l . l2009l ) , the trapping of decimeter-sized boulders in turbulen t eddies and the 



subse quent gravitational collapse of swarms of these trapped boulders (jJohansen et al 



20071 ). the trapping of particles in a pressure b ump caused by the evaporation front of 



water ( Kretke &: Lin 



2007 



Brauer et al 



2008b!) and many other scenarios. However, the 



correct modeling of any of these requires detailed knowledge of the collisional physics, 
and these models have so far relied on either simplified input phyisics or simplified initial 
conditions. 



Because of their complexity, collisional evolution models have to make simplifying as- 
sumptions about the outcome of dust aggregate collisions, for example that collisions 
always result in sticking, or otherwise use simple recipes for the collisional outcome. 
Ideally, one requires detailed knowledge of the outcome of every collision. But modeling 
this microphysics within an evolution model is simply impractica l. There are com- 
puter programs that model these individual collisions in detail (e.g.. iDominik &: Tielens 
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(|l997l ). ISuvama et alj (|2008l ): Geretshauser et al, in prep.), but each colhsion model 
takes anywhere between hours and weeks to complete on a computer. They are there- 
fore impractical to use at run-time in a model that computes the overall time-dependent 
evolution of the dust size distribution inside protoplanetary disks. These collision models 
themselves often depend on poorly known input physics. 

Another approach to obtaining the collisional outcome of dust aggregates is to model 
these collisions in the laboratory. From the many experiments that have been performed, 
a picture emerges of the outcome of dust aggregate collision under a variety of conditions 
in the protoplanetary disk (PPD). In the previous chapter, we collected data from over 
19 experiments, and compiled a set of formulae to describe reasonably well the outcomes 
of these collisions in such a way that they can be used as input to models that address 
the temporal evolution in the dust size distribution. 

In this chapter, we directly rely on the outcome of these laboratory experiments to model 
the dust aggregate size distribution. As described in Chapter [31 we have produced a 
mapping of all available collision experiments for silicate-like particles. The velocity 
range of these experiments is also sufficiently wide to cover various disk models that 
roughly correspond to the conditions at 1 AU in the PPD. For details of the collisional 
mapping, we refer to Chapter [3] but summarize the elements of our new collision model 
in Sect. 3.1. 

We build this collision kernel into a Monte Carlo code for modeling the size and poros- 
ity distribution of dust in a protoplanetary disk (jZsom fc Dullemondl (j2008l ). hereafter 
ZsD08). The outcome of our laboratory-driven dust coagulation model is difficult to pre- 
dict a priori since the key variables involved depend on a non-trivial interplay between 
the collision kernel (Chapter [3]) and the velocity field. We can, however, propose two sce- 
narios. In the first, particle growth proceeds beyond the meter-size barrier, all the way 
to forming planetesimals. In the second scenario, growth terminates at an intermediate 
size. In this case, additional growth to planetesimal sizes may proceed by the means 
of th e concentration and^ subsequent gravitational collapse of these particles f Johansen 
et al., 



2007 



Cuzzi et al.l. 



200. 



Thus, our model providea the starting conditions for 
these concentration models. We emphasize, however, that in this work we do not in any 
way 'optimize' the outcome by laboriously scanning all the parameter space or treating 
environments that may be more co nducive to growth, such as nebu l a pre ssure bumps 



2007 



Lvra et al.l . \2Q0% . These are 



or the trapping of dust in vortices (iKretke Lin , 
obvious extensions of this present work. But by considering the sensitivity of a few key 
parameters (e.g., gas density, and turbulence strength) to the outcome of the growth °o 
process, we obtain a picture of where the arrow of coagulation typically points to in ° 
protoplanetary environments: pebbles, boulders, or planetesimals. 
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In this chapter, we describe the three nebulae models used in this work and the sources 
of relative velocity between the aggregates (Sect. I4.2p . how we developed our coagula- 
tion/fragmentation model of Chapter [3] into a Monte Carlo code (Sect. 14. 3p . and our 
first results (Sect. 14. 4p . We also test the sensitivity of these results to variations in 
gas density, the velocity field, and other key model parameters. Section [4.51 reflects the 
importance of our result to planetesimal formation and provides suggestions for future 
experiments. Finally, Sect. 14.61 lists our main conclusions. 



4.2 The nebulae model 
4.2.1 Disk models 

We briefly describe the disk models considered in this chapter. 



The low density model: Resolved millimeter er nission maps of protop l aneta ry disks 
seem to indicate a shallow surface density proflle (lAndrews Williama (j2007l )) given 
by Sg(r) oc r~^'^. The systematic effects of some of their assumptions, such as the disk 
inclinations or the simplifl ed treatment of the t emperature distribution, may produce 
steeper proflles. Therefore, iBrauer et al.l (l2008al ) adopted the following profile: 



45 



g (_!_\ 
cm2 VAU/ 



-0.8 



(4.1) 



Here we assumed that the central star is of solar mass, the disk extends from 0.03 AU 
to 150 AU, and the total mass of the disk is 0.01 Mq. Assuming that the pressure 
scale-height is Hp = 0.05 x r and the vertical structure is Gaussian, we obtain: 



2TTHr, 



(4.2) 



the density at 1 AU in the midplane {z = 0) being 2.4 x 10 g cm ^, approximately 
two orders of magnitude lower than the minimum mass solar nebulae (MMSN) value. 



M MSN model: The n iinim um mass solar nebulae model (MMSN) was introduced 
by IWeidenschillinel (|l977bl ) and iHavashi et al.l (jl985l ) . Prom the present state of the 
Solar System today, it is possible to infer a lower limit to the mass in the solar nebulae 
from which the planets were formed. The model assumes that the planets were formed 
where they are currently located (no migration included). It also assumes that all the 
solid material presented in the solar nebula had been incorporated into the planets. The 
loss of solid material because of radial drift is not taken into account. Despite these 
uncertainties, the MMSN model is frequently used as a benchmark. The surface density 
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of the MMSN disk is given by 

which corresponds to a total disk mass of 0.01 Mq between 0.4 and 30 AU (between the 
orbits of Mercury and Neptune). Assuming that the vertical structure of the gas follows 
a Gaussian distribution, we infer a midplane density at 1 AU of 1.4 x 10~^ g cm~^. 



The high density model: DeschI (|2007l ) introduced a 'revised MMSN model' by 
adopting the starting positions of the planets in the Nice model of planetary dynamics 
(jTsiganis et al.l . l2005l ) thus taking into account planetary migration. The revised MMSN 
model predicts that the Solar System was initially in a far more compact configuration 
and its surface density profile is given by 



Sg(r) = 5.1 X 10^ 



g (_!_\ 
cm2 VAU/ 



-2.2 



(4.4) 



This model is consistent with that of a d eer et ion disk that is being photoevaporated by 
the central star. Although the model of iDeschI ()2007l ) was defined for the outer Solar 
System, we extrapolate the profile to 1 AU to cover a broad range of surface density 
values in our calculations. Assuming, as in the MMSN model, a Gaussian vertical dis- 
tribution, the density at 1 AU in the midplane is 2.7 x lO"*^ g cm~'^. 



For simplicity, we adopt a midplane temperature of 200 K (isothermal sound speed of 
Cs = 8.5 X 10^ cm s~^) in all three models. 



4.2.2 Relative velocities 

We consider three contributors to the relative velocities between dust aggregates: Brow- 
nian motion, radial drift and turbulence. In the following, we discuss each of these 
sources. 

The average relative velocity of two particles with mass mi and m2 in a region of a disk 
with temperature T due to Brownian motion is 



AvB{mi,m2) 



'8kT{mi + m2) 
iTmim2 



(4.5) 



For micron-sized particles, the relative velocity is of the order of 0.1 cm s~^, but for cm- 
sized particles this value drops several orders of magnitude. Therefore, Brownian motion 
is only effective for collisions between small particles during the initial stages of growth. 
Coagulation caused by Brownian motion results in fluffy aggregates of fractal dimensions 



of around 2 and 3 (Blum et al. 



19961 : iKempf et al 



199S 



Blum et al. 



200C; Krause 
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Figure 4.1: The combined relative velocities caused by Brownian motion, radial drift, 
and turbulence for fluffy particles (vP = 20) in the three disk models for equal-sized 
particles (a) and for different-sized par ticles with a mass ra tio of 100 (b). The solid 
li ne indicates the low densitv model of iBrauer et al. I (|2008ah . Physical parameters of 
tfie disk: the distance from the central star is 1 AU, te mperature is 200 K, the density 
' the gas is 2.4 x 10~^^ g cm~^, and the turbulence parameter, a = 10"'*. The dotted 
represents the MMSN model. The density is 1.4 x 10~^ g cm^^, and the other 
parameters are the same. The dashed line corresponds to the high density disk. The 



gas density is 2.7 x 10 ^ g cm ^. 



Blum. |2004| ). In practice, no growth is caused by Brownian motion for aggregates larger 
than 100 micron. 

The second contributor to the relative velocity is turbulence. The relative velocity of 
aggregat es produced by the random motion of tu rbulen t eddies was calculated n umer- 



Volketal 



(|l980h 



Mizuno et al. 



(Il988h. and 



Markiewicz et al. 



(|l99lh . We 



ically by 

use the closed form expressions pr esented by lOrmel &: Cuzzi (|2007l ). We assume that 
turbulence is parameterized by the Shakura &: Sunyaev (jl973 ) a parameter 



l^T = OtCgHq, 



(4.6) 



where vt is the turbulent viscosity, Cg is the isothermal sound speed, and Hg is the 
pressure scale height of the disk. The value of the a parameter reflects the strength 
of the turbulence in the disk. Typical values of a in this chapter range between 10~^ 
and 10^^. The turbulent relative velocity is a function of the stopping times of the two 
colliding particles. The stopping time (or friction time) is the time the particle needs to 
react to the changes in the motion of the surrounding gas. As long as the radius of the 
particle is smaller than the mean free path of the gas (a < fAmfp), the particle is in the 
Epstein regime, where the stopping time is (jEpsteinl (jl924l )): 



3m 



'iVthPgA 



(4.7) 



where m and A are the mass and the cross-section of the particle, and pg and fth are the 
gas density and the thermal velocity. At high gas densities where the mean free path 
is low or in the case of larger particles, the first Stokes regime applies and the stopping 
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Figure 4.2: The Stokes number as a function of the particle radius in the three models. 
The parameters of the dust for all of the models are the following: monomer radius is 
ao = 0.75 /im, material density is po — 2 g cm~^, and = 1. 



time is 



tst 



3m 



(4.8) 



4VthPgA 9 Amfp 

In the first Stokes regime, the stopping time is independent of the gas particle relative 
velocity as well as the gas density. This regime can be used as long as the particle 
Re ynolds number is smalle r than unity. The particle Reynolds number is calculated to 
be (|Weidenschilling| (|l977al )) 



Rer. 



2aAv, 



pg 



(4.9) 



where At;pg is the relative velocity between the particle and the gas, and t] is the gas 
viscosity. For particles outside the Epstein regime, we can assume that the system- 
atic velocity (radial drift) dominates over the random velocities (turbulence); therefore, 
Afpg ^ V£), where v_d is the drift velocity of the particle, defined in the next paragraph. 
The particle Reynolds number never exceeds unity in our simulations. Therefore, we do 
not include additional Stokes regimes. 

Radial drift also contributes to relative velocities between aggregates. Radial drift {vd) 
has two sources: the drift of individual particles {vd) and drift caused by accretion 
processes of the gas (vda), thus the total radial drift velocity i s vn = Vd + Vdn- 



The radial drift of individual dust aggregates with mass m is (IWeidenschillingi (Il977al )) ^j.^ygp ssom mjou 

"00000 

2VN 



Vd 



St + l/St' 



(4.10) 5 



where St is the Stokes number of the aggregate {St = t.^O, where O is the orbital 



frequency) and vn is the maximum radial drift velocity (jWhippld (jl972l )). 
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The second contribution to the radial velocity is produ ced by the ac c retion of the gas. 
This part of the radial velocity is calculated as follows (|Kornet et al.1 (|200lh 



Vda 



(4.11) 



where Vg^s is the accretion velocity of the gas ([Takeuchi &: LinI (l2002l )). 



The relative velocity due to radial drift is then simply the difference between the radial 
velocity of particles 1 and 2. However, as the Stokes number of the aggregates is always 
smaller than 10~^ (see Sect. 14. 4p . the second term of the radial velocity (vda) can be 
safely neglected 

AVD = \VD1 -VD2\^ \vdl -Vd2\- (4.12) 



This study uses two quantities to describe the porosity of the aggregates. The volume 
filling factor is 

(l) = V*/Vtot = {A*/Af/\ (4.13) 

where V* is the volume occupied by the monomers and Vtot is the total volume of 
the aggregate, including pores, and A and A* are the surface area equivalents of these 
quantities. In this way, the filling factors also enters the definition of the friction time 
(Eqs. 14.71 and 14. 8p . The density of aggregates then follows as p = pocp, where po = 2 
g cm~^ is the material density of the silicate. In this study, we also use the reciprocal 
parameter of the filling factor, which is denoted by the enlargement parameter, ^ = (p^^. 



We illustrate the relative velocity between equal-sized and different-sized aggregates with 
^ = 20 {(j) = 0.05) in Fig. 14.11 for the disk models considered in this work. Adopting a 
threshold (fragmentation) velocity of 1 m s~^, the maximum particle size that can be 
reached in the models are 0.025 cm in the low density model, 1.4 cm in the MMSN model, 
and 1.7 cm in the Desch model. The Stokes numbers of these particles are identical in all 
three models, 4.7 x 10^^. The constant fragmentation velocity of 1 m s~^ is the typical 
velocity at which silicate particles fragment. In our collision model this is not the case 
for all combinations of mass ratio and porosity (Chapter [3]), but the m s~^ threshold 
remains a useful proxy for the point where fragmentation processes become important. 

Figure 14.21 shows the Stokes number as a function of particle radii in the three models. 
Initially, particles are in the Epstein regime, where the stopping time, thus the Stokes 
number, depends on the gas density. When the particles enter the Stokes regime, the 
stopping time becomes independent of the gas density (see Eq. 14. 8|) . One can see that 
particles in the Desch model are in the Stokes regime at a Stokes number of 4.7 x 10~^ 
(when the particles have relative velocities of 1 m s"'^), while the aggregates in the 
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MMSN model are close to it, which explains why the maximum particle size is almost 
the same in these two models. 



As discussed in 



Qrmel &: Cuzzil (|2007l ). particles are initially in the 'tightly coupled 



particle' regime, where the eddies are all of class I type, meanin g that the t urnov er 



time of all eddies is longer than the friction time of the particles (IVolk et al 



Upon entering a class I eddy, a particle therefore forgets its initial motion and aligns 
itself with the gas motions of the eddy before the eddy decays or the particle leaves it. 
This regime is presented in Figs. 14.1b and b. Different-sized particles are found in this 
relative velocity regime as long as their masses are less than 10~® g in the low density 
model, 10^^ g in the MMSN model, and 10~^ g in the high density model assuming fluffy 
particles = 20). If the particles leave this regime and enter the 'intermediate particle' 
regime, their relative velocity increases. This transition affects the particle evolution, as 
discussed in e.g., Sect 14.4.31 



4.3 Collision model and implementation 



We use a statistical or 'particle in a box' method to compute the collisional evolution, 
that is, we assume that all particles are homogeneously distributed within a certain vol- 
ume (the simulation volume). In reality, the particles could however leave the simulated 
volume or new particles could enter from outside because of radial drift or random mo- 
tions (turbulence and Brownian motion) . Since we do not resolve the spatial dependence 
of the aggregates, we simply assume that local conditions hold during the run. The gas 
and dust densities are kept constant and particles cannot leave or enter the simulation 
volume (hereafter 'local approach'). 



4.3.1 Short overview of the colhsion model 



Many laboratory exp eriments on dust aggregate collisions have been performed (see 
Blum &: Wurrnl (120081) ). The growth begins as fractal growth and we use the recipe of 



Qrmel et al 



(|2007l ) to describe this initial stage. However, once aggregates have re- 
structured into non-fractal, macroscopic aggregates (e.g., > 100 //m in size), laboratory 
experiments show that the collisional outcomes become very diverse. In this regime, 
many new experiments were performed with dust aggrega tes consisting of 1.5 mn di- 
ameter Si02 monomers of either high porosity (j) = 0.15 (iBlum &: Schrapleii (|2004l )). 
or intermediate porosity (0 = 0.35). Chapter [3] compiled 19 experiments with different 
aggregate masses, collision velocities, and aggregate porosities. 
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Figure 4.3: The colhsion types considered in this chapter. We distinguish between 
similar-sized and different-sized particles. Some of the collision types only happens for 
one of the mass ratios. Grey color indicates that during the given collision type the 
particle is compact or part of the mass is compacted. 



From these experiments, we identified nine different collisional outcomes involving stick- 
ing, bouncing, or fragmentation (see Fig. 14. 3p . The occurrence of these regimes de- 
pends mainly on aggregate masses and collision velocities. However, it also depends on 
both the porosity of the particles and the critical mass ratio. For example. Chapter [3] 
finds fragmentation in collisions between a porous aggregate and a solid wall, whereas 
Langkowski et al. (2008) find sticking of a porous projectile by penetrating a target that 
is also porous. Likewise, Heii^elmann et al., in prep, detect the bouncing of two similar- 
sized, porous dust aggregates, while Langkowski et al. (2008) uncover sticking of the 
same velocity where one collision partner (target) is significantly bigger. To address the 
importance of the mass ratio and porosity, we identified eight different collision regimes 
(look-up tables) based on a binary treatment of porosity and mass ratio i.e., (i) similar- 
sized or different-sized collision partners and (ii) porous or compact collision partners. 
The additional distinction between target, which we always define as the heavier col- 
lision partner, and projectile then defines eight different collision regimes. We denote 
these regimes as 'pP' (porous projectile, porous target; target significantly bigger than 
the projectile) and 'pc' (porous projectile, compact target; target of similar size than 
the projectile). 

In Chapter [3l we classified each of these 19 experiments into one or more of these eight 
regimes (see Fig. 10 of Chapter [3]). Based on extrapolation of experimental findings, 
we decide in which mass and velocity range collisions result in sticking, bouncing, or 
fragmentation. These results are presented in Fig. 11 in Chapter [31 

It should be noted that the critical mass ratio of the equal-size (e.g., 'pp', 'cc') to the 
different-size regimes (e.g., 'pP', 'cC') is ill-constrained by experiments. Therefore, we 
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use critical mass ratios of = 10, 100, and 1000 to explore the effect of this parameter. 



4.3.2 Porosity 

Chapter [3] defined a binary representati on of the porosi t y, par ticles are either porous or 



compact. Following the simple model of IWeidling et al.l (120091 ) . we include a continuous 



transition between these two 'phases'. These authors showed that the compaction of 
particles caused by bouncing can be described by porous and compacted sites on the 
surface of the aggregate. A site of the aggregate is porous if it has not yet experienced 
any collisions (e.g., bouncing); a compacted site has encountered at least one collision 
but any additional collision at that part of the surface cannot change the porosity of 
the site. We define the probability of hitting a passive site of the aggregate with the 
experssion 

(pc - <Pp 

where (p is the volume-filling factor of the aggregate, (pc is the critical porosity {(pc = 0.4, 
see Chapter [3]), and (pp is the volume filling factor of the porous site, which is chosen to 
be 0.15. If (/) is between 0.15 and 0.4, a random number decides whether the particle 
collided with a porous or a compact site. This treatment of the porosity ensures a con- 
tinuous transition from porous to compact aggregates. 



During the initia l hit fc stick (SI) phase, particles are in t he fractal growth regime 



Ossenkopj (jl993l ). iBlum et al.l (l2000l ). iKrause fc BlumI (120041 )). Particles grow initially 
because of Brownian motion and later due to turbulence. The structure of the aggregate 
depends on whether the collision happened between a cluster and a monomer (PCA) or 
between two clusters (CCA). The latter type of collision produces fluffy aggregates with 
fractal dimensions of 2, while the former leads to mor e comp act structures of fractal 
dimension 3. The hit and stick recipe of lOrmel et al.l (|2007l ) attempts to interpolate 
between these two fractal models. At one point, collisional energies become high enough 
to invalidate this assumption. This o ccurs when the collis i onal energy is five times high er 



Bhim fc Wurm (1200' 



than the rolling energy of monomers (iDominik Sz Tielend (|l997l ) , 
The internal structure then becomes homogenous. In our model we assume that once 
the fractal growth caused by SI is over, bouncing with compaction (Bl) restructures 
the aggregates producing compact structures of fractal dimension 3. Since our model 
cannot follow the exact shape of the particles, we assume that the aggregates are spheres 
and can be describ ed with a single density (p) thus neglecting the effects of e.g., the 
'toothing radius' of lOssenkopa (jl993l ) or the craters forming during penetration (S3). 
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4.3.3 The Monte Carlo method 



Using the expressions for the relative vefocity, the colhsional cross-section between the 
dust particles, and the colhsional outcome, we solve for the temporal evolution of the 
dust size distribution. Traditionally, the Smoluch owski equation is used to determ ine the 
evolution of the mass distribution function fe.g . iDullemond &: Doininild (|2004^■ Dulle- 



Brauer et al. 



(j2008al )). The continuous 



mond & Dominik (|2005l ^. kanaka et alJ hom . 
form of the Smoluch owski equation used in these works lacks the stochasticity of the 
coagulation problem (jSafronovl (Il969l )). All bodies of mass m will grow in the same way 
thus the spatial and temporal fluctuations of the particle ensemble are averaged out. In 
reality, however, particles with similar masses might follow a different evolutionary path 
depending on which other particles they collide with. The collision model typically used 
in these works is, by necessity, rather simple because in the Smoluchowski formulation 
the collision and time evolution steps are linked together. These collision models consist 
of sticking and fragmentation and only the mass of the particles is followed. The advan- 
tage of such a model is that it is not computationally too expensive: the entire disk can 

be practically modeled. 
I I 



Ormel et al. (|2007l ) introduced a new Monte Carlo method to solve for the mass and the 
porosity distribution function simultaneously. Their collision model consists of sticking 
and compaction; ZsD08 also added a simple fragmentation model. Although these mod- 
els are more detailed, one can see that they still lack the full complexity observed in 
"the zoo" of laboratory collision experiments. 

The MC-approach used in this study was previously presented by ZsD08. It can be 
characterized by two key properties: (1) the number of MC-particles (also referred to as 
representative particles) is kept constant; (2) the method follows the mass of the particle 
distribution. 



Property (1) is required to ensure statistics. Because of the v-^V noise o f MC-methods, a 
large f luctuation in N would severely affect the accuracy of the method (|Ormel &: Spaans 
(|2008l )). The second property states that our primary interest lies in the particles that 
contain most of the mass of the system. Moreover, it has been shown that following the 
particle's mass distribution - rather than its number distribution - is also a prerequisite 
to preser ving a close correspondence with systems that experience strong growth f Ormel 



& Spaans ([2003)). 



Property (2) ensures that the MC method samples only the parameter space where a 
signicant portion of the total dust mass is. However, this is not always desirable. For 
instance, radiative transfer calculations require as input the surface area distribution of 
the aggregates, which determines the opacity. If most of the particle mass is contained 
in large particles (which are not observable), the number of small particles (which could 
contain most of the surface area and determines the IR appearance of the disk) might 
be resolved with poor statistics. But if we are interested in following the evolution of 
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the dominant portion of the dust, then MC methods naturaUy focus on these parts of 
phase space. 

A necessary condition for the ZsD08 method to work is that the number of the repre- 
sentative particles is much less than the number of actual aggregates present in the 
system being considered - a condition that is safely met in any of our simulation runs. A 
representative particle will then collide only with the non-representative particles, whose 
distribution is assumed to be identical to that of the representative particles. We refer 
to ZsD08 for details of the precise implementation and accuracy of the method; here we 
concentrate further on how the method operates in the new collisional setup. 



The collision kernel is defined as the product of the cross-section of the colliding particles 
and their relative velocity: 

Ki^k = (Ti^k^ViM, (4.15) 

where the index i corresponds to the representative particle and k is the index of the 
non-representative particle. The kernel is proportional to the probability of a collision. 
The value of ATj^fc is calculated for every possible particle pair, and random numbers 
determine which of the collisions occur first and at which time interval. 



The above properties and conditions specify the essence of the ZsD08 method: one 
of the two collision particles is a representative particle and, because of property (1), 
only one of the collisional products becomes the new representative particle. Because 
of property (2), the choice for the new representative particle is weighed by the mass of 
the collision products. A very helpful analogy here is that of the representative 'atom', 
which is contained within the representative particle. The choice of new representative 
particle after the collision is then proportional to the probability of the representative 
'atom' ending up in the collision products. If, for instance, a collision leads to the 
production of an entire distribution of debris particles, the probability that a particular 
debris fragment becomes the new representative particle is proportional to the likelihood 
of this fragment containing the representative 'atom'. 



4.3.4 Implementation of the collision types 



We describe the implementation of the collision model using the representative 'atom' , 

^ Of A^jSUSp SSDLU -LUJOU 

concept. We refer to Chapter [3] for details of the various collision types mentioned below, o b o o o 



p IT) ^ G 

d d d 



Hit & stick (SI), sticking through surface effects (S2), penetration (S3): 

All three of these collision types cause sticking and increase the mass of the aggregate 
by that of the projectile, but the porosity changes in a different manner (see Chapter 
(31) . The new mass of the representative particle i is then the sum of the original particle 
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masses, mj^new = n^-i + mk, where rui is the mass of the representative particle and 
is the mass of the non-representative particle. 



Mass transfer (S4): In the case of mass transfer (S4), a certain percentage of 
the mass of the projectile sticks to the target, while the leftover mass of the projectile 
fragments according to a power-law distribution (see Chapter [3]) . 

There are two situations to consider: 



1. The representative 'atom' is part of the target. The mass of the new aggregate 
will be the mass of the original aggregate plus the transferred mass from the non- 
representative particle (mj^ncw = "mi + "T-transi where mtrans is the transferred mass 
calculated according to Chapter [3]) . 

2. The representative 'atom' is part of the projectile. We again have two situations. 

(a) The representative 'atom' will be transferred to the non-representative parti- 
cle. The mass of the new representative particle will be the mass of the non- 
representative particle plus the transferred material (mj^new = w-fc + "itrans)- 
The probability of transferring (removing) the representative atom from the 
projectile is simply P = Jn-trans/^T-j, the ratio of the transferred mass to the 
mass of the projectile. 

(b) The representative 'atom' remains in one of the fragments. The probability 
of this event is P = (m^ — mtrans)/™'^) the ratio of the fragmented mass 
to the original mass of the representative particle. As discussed in Chapter 
[31 the fragments follow a power-law mass distribution. The distribution is 
defined by the maximum mass of the fragments, which is a function of the 
relative velocity and the total mass of the fragments. The total mass of the 
fragments is mi — mtrans- We randomly choose from the fragment distribution 
to determine the new mass of the representative particle (to find which of the 
fragments will contain the representative 'atom'). 



Bouncing with compaction (Bl): In this process, particles collide and bounce. 
Bouncing itself does not change the mass of the particles, b ut it compactifies the m 
according to Chapter [3l As observed in laboratory experiments (jWeidling et al.l (j2009l )). 
there is a small probability (Pfrag = 10~^) that the bouncing particle will break apart. 
If this happens, we break the particle into two equal-mass pieces. 



Bouncing with mass transfer (B2): From the implementation point of view, this 
is similar to mass transfer (S4). The recipe to define the new representative particle is 
as in mass transfer (S4). The difference is that the projectile does not fragment during 
the collision, and that the porosity changes differently (see Chapter [3|). 
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Pragmentation (Fl): Fragmentation only happens between similar-sized aggregates 
in the 'pp ' and 'cc ' regimes. The fragments follow a power-law mass distribution, where 
the maximum mass of the fragments is determined by the relative velocity of the particles 
and the total mass that goes into the fragments (Chapter [3]). We randomly choose from 
this distribution to determine the new mass of the representative particle. 

Erosion (F2): Erosion (F2) happens between only different-sized particles. During 
the collision, the projectile "kicks out" pieces from the target aggregate. These pieces 
follow a power-law distribution (see Chapter [3]). We have to consider two cases: 

1. The representative 'atom' is in the target. Again, we have two possibilities. 

(a) The representative 'atom' remains in the target after the collision. The mass 
of the new particle will be ?n.j,new = — w-er, where ruer is the eroded mass. 
The probability of this event is P = (mj — mer)/?7ii, which is the ratio of the 
left-over mass (which does not erode) to the mass of the original particle. 

(b) The representative 'atom' originates in the eroded particles. Since the eroded 
particles follow a power-law distribution, we randomly draw from this distri- 
bution to determine the new mass of the representative particle. The likeli- 
hood of this event is the ratio of the eroded mass to the original mass of the 
particle {P = mgr/mj). 

2. The representative 'atom' is part of the small particle that caused the erosion. As 
the particles do not stick and the small particle does not fragment, the represen- 
tative particle remains unaffected. 



Fragmentation with mass transfer (F3): The porous particle becomes destroyed 
by the compact one and transfers a certain amount of mass to the compact particle, 
fragmentation with mass transfer (F3) only happens in the 'cp' regime. We again have 
two possibilities. 

1. The representative 'atom' is part of the compact particle. In this case, the rep- 
resentative 'atom' cannot leave the particle. The new mass of the representative 
particle will be mj^new = "ij + mtrans, the sum of the original mass plus the trans- 
ferred mass. 

2. The representative 'atom' was part of the porous aggregate. 

(a) The representative 'atom' is part of the material that is transferred to the 
compact particle. In this case, the new mass of the particle will be that 
of the compact (non-representative) particle plus the transferred material 
("ii,new = + "^trans)- The probability of this event is P = mtrans/"^i- 
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(b) The representative 'atom' is among the fragments. As before, the mass dis- 
tribution follows a power-law and we draw randomly from this distribution 
to determine the new mass of the representative particle. The probability of 
this event is P = (mj — mtrans)/'"^i- 

4.3.5 Evolving the particle properties in time 

We summarize how the particle properties evolve in time using the aforementioned 
kernel. We begin with the size and porosity distribution of the particles at a given 
time, t. At t = 0, we provide the initial size and porosity distribution (see Sect. I4.4.ip . 
Knowing these: 

• We calculate both the cross-sections of all possible collision partners and their rel- 
ative velocities using the equations described in Sect. 14.2.21 Both sets of quantities 
are used to determine the collision rates between the particle pairs. 

• By using random numbers, we identify from the collision rates the representative 
particle involved in the collision, the non-representative particle it collides with, 
and the time at which the collision takes place {t + At). 

• Knowing the masses (mass ratio) and porosities of the collision partners, we iden- 
tify in which of the eight regimes the collision takes place (e.g., 'pP', or 'pC, 
etc.). 

• Next, we identify which of the nine collision types materializes (Fig. 14. 3p using the 
relative velocity of the particles and the mass of the projectile (see Chapter E]). 

• Based on the collision recipe described in Chapter [3] and Sect. I4.3.4| the new mass 
and new porosity of the representative particle is calculated and the new size and 
porosity distribution of the particles at time t + At is obtained. 

• In the final step, we update the collision rates. 

4.3.6 Numerical issues 

As mentioned in ZsD08, a sufficiently high number of representative particles is needed 
to properly reproduce the physics of the collision kernel. We performed simulations with 
an increasing number of representative particles and found that for more than 200-300 
particles the results of the simulations do not change significantly. For all simulations 
described in the following sections, 500 representative particles are used and we average 
the results of 20 simulations to decrease the numerical uncertainty of the code. 
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The required computational time strongly depends on the collision rate of the particles 
thus determined by the dust density and the relative velocity (the a turbulence param- 
eter mostly) and the length of the simulation. On a 2.83 GHz CPU, performing the 
simulations on a single core, the CPU time varies between nine hours (the low density 
model with a = 10~^) and three days (the high density model with a = 10~^). Both 
simulations modeled 10^ years of particle evolution. The high density simulation with 
a = 10~^, critical mass ratio of 100, and t = 10'' years of evolution takes twelve days to 
complete. 



4.4 Results 

4.4.1 Initial conditions and setup of simulations 

All simulations begin with silicate monomers of 1.5 //m diameter and 2 g cm~^ material 
density (monodisperse size distribution) . We simulate the dust evolution at the midplane 
of our disk models at a distance of 1 AU from the central star. The gas density is 
obtained from the disk models described in Sect. 14.2.11 We assume a typical 1:100 
dust-to-gas ratio. We follow the history of each collision: the mass and porosity of the 
colliding particles, their relative velocity, the occurred collision type, and the new mass 
and porosity of the particles. In this way we are able to reconstruct the history of the 
dust evolution. 

The parameters that we vary in this study are the gas density pg and the turbulence 
parameter a. We also treat the critical mass ratio as a free parameter to explore its 
effect on the dust evolution. 

We provide a detailed description of the low density model with a = 10~^ and critical 
mass ratio of 100 in Sect. 14.4.21 We then compare this with the MMSN model and 
the high density model using the same turbulence parameter and the critical mass ratio 
(Sects. 14.4.31 and I4.4.4p . In Sects. 14.4.51 and 14.4.61 we discuss the effects of changing the 
turbulence parameter and critical mass ratio by comparing those results with the two 
example runs. 



4.4.2 The low density model 
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In this disk model, the gas density at 1 AU is 2.4 x 10 g cm ^, 



the turbulence 2 ° ° ° ° 
parameter is a = 10~^, and the critical mass ratio is r„i = 100. As shown in Fig. 5 
14. H the particles reach the fragmentation velocity (1 m s^^) at sizes smaller than a 
millimeter because the particles in low gas density environment decouple from the gas 
at these small radii. 
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Figure 4.4: (a) The evolution of the mass distribution, (b) enlargement parameter 
distribution, and (c) the collision frequency of the nine different collision types in the 
low density model with a = 10~^ and critical-mass ratio of 100. The x-axis shows 
the time. The y-axis of the (a) and (b) figures show the logarithmic mass and the 
linear enlargement parameter, respectively. The contours represent the normalized mass 
density and the mass weighted enlargement parameter. The black lines represents the 
average of the mass and enlargement parameter at a given time. The y-axis on the (c) 
figure represents the nine collision types. Each stripe shows the total collision rate of 
the collision types. Two distinct phases can be distinguished. During the initial 300 yr, 
particles grow by hit & stick (SI), after that the evolution is governed by bouncing with 
compaction (Bl). The white lines indicate how long our 'local approach' assumption 
remains valid (discussed in Sect. 14. 3p . 
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Figure 4.5: The collision history of the eight regimes in the low density model for 
a = 10"''. The x-axis is the relative velocity, the y-axis shows the projectile mass. The 
different collision types, their border lines, as well as the areas covered with laboratory 
experiments (grey) are plotted. A relative velocity - mass grid is created and in these 
grid cells we calculate how many collisions happened until the 'local approach' assump- 
tion is valid (4 x 10^ yr). This is represented by the colors: yellow and red indicate 
a high collision frequency. The two dotted lines in the 'cc' regime are evolutionary 
tracks. Assuming a constant (40%) volume-filling factor, the relative velocity between 
equal-sized particles (left curve) and particles with a mass ratio of 100 (right curve) 
can be calculated. The collisions in the simulation should occupy the parameter space 
between these two lines. The small deviations are due to the volume filling factor not 
being exactly 40% during the simulation. 
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Figure shows the evolution of the mass distribution (a), the porosity distribution 
(b), and the cohision frequency of the various colhsion types (c). The x-axis shows the 
time on a logarithmic scale. The y-axis of Fig. 14.4b . b shows the mass and enlargement 
distributions, respectively. The intensity of the color reflects the number density of 
representative particles, which, as explained in ZsD08, indicates the mass density of the 
distribution. Thus, in Fig. 14.4b the intensity levels directly reflect the mass density, 
while in Fig. 14.4b the colors indicate the mass-weighted enlargement parameter. The 
black lines show the average of these quantities over the particle distribution. The y-axis 
in Fig. 14.4b represents the nine collision types used in this chapter. Every stripe shows 
the total collision rate of the collision types at a given time. 

Figure [^31 represents the collision history in the eight collision regimes. The x-axis is the 
velocity, the y-axis shows the mass of the projectile. A mass- velocity grid is created and 
we calculate how many collisions occurred within each grid cell during which our 'local 
approach' assumption is correct. The different collision types and their border lines, as 
well as the areas covered by the laboratory experiments (indicated by grey colors) are 
plotted. For more details on the experiments, we refer to Chapter [3l 

In the 'cc' panel, we indicate two curves with dotted lines. These curves are evolution 
tracks. The left curve is obtained by calculating the relative velocity between equal-sized 
particles with an enlargement parameter of 2.5 (volume-filling factor of 40%). The right 
curve represents the relative velocity between particles with a mass ratio of 100. These 
two curves serve as a guide to our results, as collisions should occur between these two 
curves in the 'cc' panel. The lower part of the left curve, where the relative velocity 
decreases with increasing mass, is a sign that relative velocities between equal-sized 
particles are dominated by Brownian motion. For higher masses, the relative velocity 
is dominated by turbulence. These curves do not precisely match the contours because 
we assumed a constant enlargement parameter of 40% when calculating the evolution 
tracks, whereas ^ is a free parameter in the simulation. 

4.4.2.1 Early evolution 

We discuss the evolution of the distribution functions until the 'local approach' assump- 
tion becomes invalid (4 x 10^ yr). The long-term evolution of the dust is discussed in 
Section iXMl 

We distinguish between two distinct phases. During the first 300 yr, particles grow by the 
means of the hit &: stick (SI) mechanism. The second phase is dominated by bouncing 
with compaction (Bl); the particles leave the SI regimes. During this phase, the mass 
of the particles slowly decreases and the enlargement parameter asymptotically reaches 
a minimum value of 2.23. As discussed in Chapter [3l keeping the bouncing velocity of 
a particle constant, the porosity of the aggregate will asymptotically reach a maximum 
value, (/'max (see Chapter [3]). The relative velocity of a particle is a function of the friction 
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time (Eq. 14. 7|) . which depends on the ratio of the mass to surface area, m/A. Since 
particle growth is halted at this point in the simulation (m remain constant), only a 
decrease in A caused by compaction can further increase the velocity between particles. 
The particle radius can decrease until either 0max for the given relative velocity is reached 
or particles reach the maximum compaction possible. The latter limit, random close 
packing (RCP), corresponds to an enlargement parameter of 1.6 (volume- filling factor 
of ~60%). 

We find that fragmentation does not play a role during the evolution of these particles 
indicated by Fig. 14.4b . As can be seen in Fig. 14. 5^ their evolution is halted by bounc- 
ing before the particles are able to reach the fragmentation barrier. The two dominant 
collision types are hit & stick (SI) and bouncing with compaction (Bl). 



4.4.2.2 Termination of growth 

As we can see from Fig. 14.51 sticking at higher energies than the hit & stick (SI) limit 
is possible only inside the 'pP' regime. As soon as we no longer have collisions inside 
this regime or the SI regimes, the growth is halted. There can be two reasons why this 
occurs 1.) All particles become compact, i.e., there are simply no collisions in the 'pP' 
regime. 2.) The width of the particle mass distribution is smaller than the critical mass 
ratio (r^), such that all collisions take place in the equal-size regimes (e.g., 'pp', 'pc'). 

In the case of the current simulation, the small particles have been 'consumed'. Once the 
heavy particles grow into the bouncing with compaction (Bl) area of the 'pp' regime, 
their growth in the 'pp' regime stops. The heavy particles collect the small ones via 
collisions in the 'pP' regime and by doing so, the width of the distribution is reduced 
to a value smaller than r^- Therefore, before particles are able reach the fragmentation 
barrier, growth is halted. Because of Bl, particles become compact and collisions in the 
'cc', 'cp', and 'pc' regimes appear. 



4.4.2.3 Long-term evolution 



Before discussing the long-term evolution of the distribution functions, we must consider 
how long our initial assumptions ('local approach' and constant gas density) hold true. 
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Using Eq. 14. IH we calculate that a particle with Stokes number 10~^ drifts a distance of 
1 AU in roughly 4 x 10^ yr. This is the drift timescale beyond which the 'local approach' 5 
assumption (discussed in Sect. 14. 3|) is no longer valid and particles become separated 
from each other on this timescale. 
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Another process by means of which particles separate is viscous spreading. We assume 
that the viscous timescale of the disk at 1 AU is given by: 

ivis = r^UT, (4.16) 

where r is the distance from the central star (1 AU), vt is defined in Eq. 14.61 The 
viscous timescale in our model, using a = 10~^, is of the order of 10^ yr. 

One has to consider the results of the simulation with caution for longer times than the 
drift or viscous timescales. The equilibrium or final state of the particles is reached when 
mass decrease during bouncing with compaction (Bl) and bouncing with mass transfer 
(B2), and mass increase by hit & stick (SI) and sticking through surface effects (S2) 
are in equilibrium. In other words, the final state is reached when the evolution of the 
average mass and the enlargement parameter can only be determined by the stochastic 
fluctuations of the simulation. We find that the equilibrium state of the particles is 
hardly reached within these timescales. Upon neglecting these warnings, we find that 
the equilibrium state of the dust is reached at t = 4 x 10^ yr. The equilibrium is reached 
between the bouncing collisions, resulting in breakage and hit & stick (SI) (see Fig. 
14.4b ). The equilibrium average mass and porosity of the particles are mfin = 2 x 10^^ g, 
and ^fin = 2.77. 

To be able to compare the distribution functions of different runs, we define some quan- 
tities using the mean of the distribution functions indicated with black lines in Figs. 
14.4b and b: max(m), the maximum of the mean mass; max(^), the maximum of the 
mean enlargement parameter; "^min, the minimum mean enlargement parameter when 
particles can no longer be compacted anymore; inoo the time when ^'mm is reached, 
which is when the time derivative of ^ becomes zero {d^/dt = 0); and max(S't), the 
maximum average Stokes number reached during the simulation. The values of these 
quantities are listed in Table [d] (model id 'Ltld-4ml00'). In this table. Col. 1 describes 
the model names, 'L' representing the low density model, 'M' being the MMSN model, 
'H' being the high density model, the letter 't' and the following number indicates the 
value of the turbulence parameter, and the letter 'm', and the number providing the 
used critical mass ratio values. Columns 2, 3, and 4 show the gas density, turbulence 
parameter, and the critical mass ratio respectively. Columns 5, 6, 7, 8, and 9 list the 
parameters defined to characterize the distribution functions. These are max(m) in Col. 
5, max(^') in Col. 6, ^min in Col. 7, tnoc in Col. 8, and finally max(5't) in Col. 9. 

4.4.3 The MMSN model 

In the MMSN model, the gas density at 1 AU at the midplane is 1.4 x 10~^ g cm~^, 
a = 10~^, and the critical mass ratio is 100. As shown in Fig. 14. H the particles become 
larger than in the low density model, because they are more tightly coupled to the gas 
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and the relative velocities are suppressed. As in the previous section, we first discuss 
the evolution of the distribution functions for as long as the 'local approach' assumption 
holds true (6 x 10^ yr in this model). 

Figure 14.61 again shows the time evolution of the mass (a) , the enlargement parameter 
(b), and the collision frequency (c). Figure H7fl shows the collision history. These figures 
depict a rather different evolution than that of the previous model. 

4.4.3.1 Early evolution 

During the fractal growth regime, we find that the collision rate of hit &: stick (SI) is 
much higher than in the low density model (Fig. 14.6b ). This is because of the higher 
dust densities. We can see from Fig. 14.71 for the 'cc' regime, that growth begins with 
Brownian motion because the relative velocity decreases with increasing particle mass 
for particle masses less than 10^^ g. As a result of these low velocity collisions, some 
particles reach enlargement parameter values of higher than 30 (volume-filling factors 
less than of 3.3%). At 200 yr, some particles grow above the border line of hit & stick 
(SI) and enter the area of sticking through surface effects (S2) in the 'pP' plot, and 
bouncing with compaction (Bl) in the 'pp' plot. Growth caused by both SI and S2 
continues until different-sized particles enter the transition regime in the 'pP' plot. One 
can see in Fig. 14.6b . that some particles reach 1 g in mass. However, when particle 
collisions enter the transition regime between bouncing with mass transfer (B2), and 
sticking through surface effects (S2) in the 'pP' plot, their masses are equalized because 
of the mass transfer that occurs during the B2 collisions and the collisions shift to the 
similar-sized regime (Bl). After roughly 10^ yr, we find that particles mostly bounce 
and become compact. The enlargement parameter reaches a minimum value of 1.85 
(54% volume filling factor), the mass distribution function slowly decreases because of 
the low probability of breakage. Collisions at this point occur mainly in the 'cc ' regime. 

A peculiar feature of Fig. 14.6b is a peak at t = 2.5 x 10^ yr, which is accompanied by 
a fast decrease in the enlargement parameter in Fig. 14.6b and an increased collision 
rate of sticking through surface effects (S2) and bouncing with mass transfer (B2) in 
Fig. 14.6b . At this point, the relative velocity produced by turbulence increases. As 
discussed in Sect. I4.2.2| particles leave the 'tightly coupled particle' regime and enter 
the 'intermediate particle' regime (see the relative velocity bump in Fig. 14.1b ). We 
calculate the growth timescale of the heaviest particle with mass M in the simulation % 



to be: 



o 




(4.17) 
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Fr, with moss tronsfer — F3 
Erosion - F2 
Frogmentotion — F1 
B. with moss transfer — B2 
Bouncing with compaction - B1 
Moss tronsfer - S4 
Penetrotion - S3 
St. through surfoce effects - S2 
Hit&stick - SI 

10' 




10" 

time [yeors] 



Figure 4.6: Same as Fig. 14.41 but for the MMSN modeL We magnify the spike of the 
mass distribution at ^ 2.5 x 10'^ yr in Fig. (a). Four phases can be distinguished here. 
InitiaUy (first 300 yr) particles grow purely by hit & stick (SI). After this, the growth 
slows down because Bouncing with compation (Bl) starts and all particles leave the hit 
& stick (SI) regime. Between 3 x 10"^ and 10^ yr, particles enter the transition regime 
between sticking through surface effects (S2) and bouncing with mass transfer (B2) 
on the 'pP' regime. Some particles reach masses of 1 g, but their masses are lowered 
rapidly by B2. The last phase is dominated by bouncing with compaction (Bl). The 
solid/dotted white lines indicate how long our 'local approach' assumptions remains 

valid (discussed in Sect. 14. 3|) . 
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Figure 4.7: Same as Fig. I4.5l but for the MMSN model. The particles are more tightly 
coupled to the gas due to the higher gas density. Therefore, they grow to larger sizes 

than in the low density model. 
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Figure 4.8: The dotted line and the '+' signs represent the growth timescale of 
the heaviest particle in the MMSN simulation with a = 10^'' and r„, = 100. As 
a comparison, we show the minimum growth timescale a particle can have in this 

simulation (solid line). 



which is illustrated with a dotted line in Fig. 14.81 As a comparison, we also calculate 
the minimum growth timescale that a particle can have (solid line), which 

where is the cross-section of the largest particle. Here, we assume that the 'swept up' 
particles have masses of M/lOO, therefore we use the relative velocity curve presented in 
Fig. 14.1b . dotted line. The effects of both the relative velocity bump and the increased 
growth rate can be seen at 0.1 g. 

The relative velocity 'boost' happens shortly after the particles enter the transition 
regime of S2 and B2 in the 'pP' plot. The heaviest particle, which experiences the ve- 
locity transition the earliest, acquires higher relative velocities, leading to a higher rate 
of collision with the other particles. As the particles are initially in the lower part of the 
S2-B2 transition regime (with masses of 10~^ g, see Fig. 14.6b ). the heaviest particle ex- 
periences rapid growth reaching masses of 30 g. The simulated timescale, however, does 
not reach the minimum growth timescale because of the bouncing with mass transfer 
(B2) collisions, which reduce the mass of the heaviest particle. The remainder of the 
particle population increase in mass because of B2, and the growth rate of the heaviest 
particle decreases. Eventually, the rapid growth of the heaviest particle is halted, and 
the growth timescale at m = 30 g becomes infinity. From this point on, the heaviest 
particle reduced in mass, and B2 equalizes the masses of the particles. 
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4.4.3.2 Long-term evolution 

We calculate the drift and viscous timescales to determine how long our assumptions of 
'local approach' and constant gas density remain valid. Assuming Stokes number 10~^ 
particles, we find that the drift timescale is of the order of 6 x 10^ yr, and the viscous 
timescale is 10^ yr. These timescales are indicated by solid and dotted white lines in 
Fig. IMl 

We find that the final equilibrium is reached at t = 2 x 10^ yr, which is longer than the 
drift and viscous timescales. The equilibrium is reached between the growth mechanisms 
of hit & stick (SI), sticking through surface effects (S2) and the destruction mechanisms 
of bouncing resulting in breakage and bouncing with mass transfer (B2). The final av- 
erage mass and porosity of the particles are mgn = 2 x 10^'^ g, ^^fin = 3.3. 

We conclude that the dust evolution is more complex in the MMSN model than in the 
low density model because the complex interaction of the velocity field and the collision 
kernel is apparent in this model. As in the previous model, bouncing with compaction 
(Bl) is the most frequent collision type and hit & stick (SI) determines the initial particle 
growth, but both sticking through surface effects (S2) and bouncing with mass transfer 
(B2) are of importance in this model. The final equilibrium is not reached within the 
drift and viscous timescales. 

4.4.4 The high density model 

The gas density in this model is 2.7 x 10^^ g cm^'^ at the midplane of the disk at a 
distance of 1 AU from the central star. The values of a, r^, and the dust-to-gas ratio 
are the same as in the previous models. 

Figure HTH dashed line, shows the relative velocity field of fluffy aggregates in this model. 
As already discussed in Sect. I4.2.2| the aggregates reach 1 m s~^ relative velocities at 
similar masses as the MMSN model, because of the Stokes drag. Therefore, we expect 
that the final aggregate sizes and masses will be similar to the particles produced in the 
MMSN model. 

Figure [^M shows the time evolution of the mass (a), enlargement parameter (b), and 
the collision frequency (c). Figure H.IOI illustrates the collision history. 

4.4.4.1 Early evolution 

As seen in Fig. I4.10| Brownian motion is the dominant source of relative velocity, as 
long as particles have masses of lower than 10"*^ g (that is an order of magnitude higher 
than in the MMSN model). Therefore, the enlargement parameter of the aggregates 
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Figure 4.9: Same as Fig. 14.41 but for the high density modeL As in Fig. 14.6b . we 
zoom in on the peak at the mass distribution. The sohd white hne indicate how long 
our 'local approach' assumptions remains valid at t = 10^ yr (discussed in Sect. 14. 3p . 



108 



Chapter 4. Bouncing barrier 



velocity rcm/s] 
10"^ 10^ 10 



velocity [cm/sl 




0' 

CD 

0=\ 

m 
c 



o 

0'^ 



O'f 



CP 

0^\ 

CO 

<z 



o 
o 

0' o 
O'f 



C 

0^ 

"o 
o 

0' o 
O'f 



0' o 

o'f 

0° " 



10"' 10^ 10 

velocity [cm/s] 



velocity [cm/s] 



Figure 4.10: Same as Fig. I4.5l but for the high gas density model. 
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is also higher than in the MMSN model. As the hit & stick (SI) collisions are more 
frequent than in the MMSN model because of the higher dust densities, the particles 
reach the transition regime between sticking through surface effects (S2) - bouncing 
with mass transfer (B2) earlier, at i = 200 yr. The peak in the mass distribution is not 
as pronounced as in the MMSN model. The relative velocity boost occurs at heavier 
aggregate masses (10~^ g, see Fig. 14.1b ) because of the higher gas density of the model. 
When the rapid growth of the heaviest particle begins, most of the projectiles are already 
in the transition regime. Here, the B2 collisions soon reduce the mass of the heaviest 
particle and narrow the mass distribution. 

In contrast to the MMSN model, the mass of the particles is not reduced because of 
the low probability of breakage in bouncing with compaction (Bl), but is kept nearly 
constant in time. This is the result of the increased collision rate of sticking through 
surface effects (S2). The S2 collision rate is increased because of the low velocity col- 
lisions, which occur when particles are in the tightly coupled regime and have similar 
stopping times. These S2 collisions occur in the 'pp' regime as seen in Fig. 14.101 These 
collisions cancel out the effect of breakage in Bl. 

The maximum Stokes number reached in this model is 3.6 x 10~^ (see Table model 
id 'Htld-4ml00'), which is lower than in the MMSN model. The growth in this model 
is halted by the bouncing with mass transfer (B2) collisions in the transition regime of 
the 'pP' panel. This shows us that particles cannot reach masses much higher than 1 
g independently from the exact value of the gas density, because at this point, particles 
enter the S2-B2 transition regime and the growth is halted. Increasing the gas density 
yet further would lead to even lower Stokes numbers. 



4.4.4.2 Long-term evolution 

The drift and the viscous timescales in the high density model are both 10^ yr. As 
seen in Fig. 14.9b . the particle masses do not change significantly after t = 10^ yr. The 
porosity is reduced by bouncing with compaction (Bl) and reaches a final value of 5.41 
at t = 3 X 10^ yr. 



4.4.5 Varying the turbulence parameter 



To explore the effects of turbulence, we perform two more simulations in each of the 
disk models. We keep the critical mass ratio fixed (100) and vary only the turbulence 
parameter (a) to obtain values of 10"'^, 10~^, and 10^^. The results are shown in Table 
14. H for the first nine models, and in Fig. 14.11b . 



The work of iBrauer et al.l (|2008al ) suggests that in situations where fragmentation limits 
the growth, a lower turbulence strength produces larger aggregates. This, of course. 
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Table 4.1: Overview and results of all the simulations. 
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In this table, Col. 1 describes the model names. 'L' stands for the low density model, 
'M' is the MMSN model, 'H' is the high density model, the letter 't' and the following 
number indicates the value of the turbulence parameter, the letter 'm' and the number 
shows the used critical mass ratio values. Columns 2, 3, and 4 shows the gas density, 
turbulence parameter, and the critical mass ratio respectively. Columns 5, 6, 7, 8, and 
9 list the parameters defined to characterize the distribution functions. These are the 
average maximum mass in Col. 5, the average maximum enlargement parameter in Col. 
6, the minimum enlargement parameter in Col. 7, the end of the compaction phase in 
Col. 8, and the average maximum Stokes number in Col. 9. 
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Figure 4.11: The maximum mean particle mass as a function of the turbulence pa- 
rameter (a) and critical mass ratio (b). 
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directly reflects the shift in the fragmentation threshold (1 m/s) to larger sizes when 
a is lower (Fig. I4.ip . In this study, it is fragmentation that balances the growth, 
producing a (quasi) steady-state. For the low density models, we do see a decrease in 
the final particle mass, but it is balanced by bouncing and not fragmentation. In the 
low density model, particles grow only in the hit & stick (SI) regimes. When particles 
leave these regimes, the growth stops due to bouncing. The border of the SI regime 
corresponds to the collision energy that is lower than 5 x £"1011, where £'roii is the rolling 
energy of monomers (see Chapter [3]). As the collision energy is i?coii = l/2/u(Aw)^, 
particles in strong turbulence leave the SI regimes at lower particle masses. 

On the other hand, the MMSN and high density models show that the maximum mass 
of the particles can even increase with a. The precise value of the max(?fi) is determined 
by the intensity of the peak in the mass-density plots (Sect. 14.4.3. ip and this may vary 
between simulations. In the 'Mtld-4ml00' model, we have argued that the spike is 
exceptionally pronounced because of the high probability of sticking through surface 
effects (S2) collisions at the initial part of the rapid growth. However, the main point 
is that in the MMSN/high-density simulations the maximum particle masses all end up 
around 1 g, independent of the turbulent strength because of the nature of the S2-B2 
transition, which occurs at projectile masses of 10~^ g in the 'pP' plot. As explained 
before, collisions in the 'pP' plot are the only way by which particles can grow after the 
hit & stick (SI) phase is finished. Thus, we require a broad distribution with a high 
growth rate. However, B2 collisions operates in the opposite way: they transfer mass 
from the target to the projectile, narrowing the distribution and decreasing the overall 
probability of the 'pP' process occuring. Thus, once B2 becomes effective, there is a 
shift from the 'pP' panel to the 'pp' panel. For the MMSN/high-density models this 
behavior is always present and the important quantities involved (i.e., relative probability 
of B2 over S2) scale with mass but not with velocity. The result is that the maximum 
masses that particles reach are ~ 1 g, which is rather insensitive to the strength of the 
turbulence. 

4.4.6 Varying the critical mass ratio 

We perform simulations in the disk models with a = 10~^ but a varying critical mass 
ratio. We explore how the dust distributions change upon using r^ = 10, 100, and 1000. 
Table St], lines 10 to 18, shows the parameters that describe the distribution functions, 
and Fig. 14.11b illustrates the maximum particle mass as a function of the critical mass 
ratio. 

By examining Table 14. H we see that using = 10 in the low density model ('Ltld- 
4ml0') results in heavier and more compact particles. The low critical mass ratio means 
that the largest particles in the different-sized regimes can sweep up the projectiles and 
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Figure 4.12: The colhsion frequencies of the 9 coUision types in the MMSN model 

with a — 10^^ and rm = 100. 



grow to larger sizes, eventually reaching the fragmentation line, where growth stops. As 
discussed in Sect. 14. 2. 2^ assuming a fragmentation velocity of 1 m s~^, the maximum 
Stokes number of the aggregates is 4.7 x 10~^, a value almost reached in this model. 

We find that there is no significant difference between the r^ = 100 and 1000 simulations 
in the low density model. The explanation for this can be found by examining the width 
of the mass distribution in the hit & stick (SI) phase. This initial phase occurs in the 
same way independently of the critical mass ratio. If the critical mass ratio is equal 
to or larger than the width of the distribution function, collisions between different- 
size particles in the 'pP' regime are inhibited. After the SI phase, the width of the 
distribution in the low density regime is approximately 100. Therefore, we do not see 
any difference when the mass threshold is shifted from = 100 to = 1000; in both 
cases, collisions occur between equal-size particles only, and these are either SI or (when 
this stage is over) Bl. 

For the high density models (MMSN/Desch), we find that the outcome is again similar: 
growth halts at ~ 0.1 g (within a factor of 10) and no clear dependence on is seen. 
For the high mass ratios, growth is always in the similar-size regime. Here, it is the gas 
density that determines the velocity, i.e., whether we have a sticking (SI) or a bouncing 
(Bl) collision. Therefore, if = 1000, the high density model produces heavier particles 
than the MMSN model (see Fig. 14.11b ). For lower r^, it is again the nature of the S2-B2 
transition regime that limits the maximum mass. 

Thus, the critical mass ratio is an important parameter because it determines the relative 
likelihood of collisions occurring in the different-size regime, which are in general more 
conducive to growth. In contrast, for simulations where B2 collisions are important 
- which have the effect of narrowing the distribution - the width of the distribution 
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will correspond to the value of the parameter, although we have also seen that the 
absolute size/mass is rather insensitive to the parameter. Overall, these arguments 
indicate that a good knowledge of this parameter is important. 

4.5 Discussion 

We have performed simulations of varying turbulence parameter and critical mass ratio 
values in three disk models with low, intermediate and high gas densities. We have 
found that hit & stick (SI) and bouncing with compaction (Bl) are the most dominant 
collision types. All simulations show evidence of long-lived, quasi-steady states. Frag- 
mentation is rarely present, but even then, for only a limited time period. The absence 
of fragmentation is caused by the bouncing collisions. 

4.5.1 The sensitivity of the results 

As presented in Sect. 14.41 the outcome of our simulations is determined by the collision 
kernel, and the relative velocity field. A significant change in one, or both can alter the 
evolution of the aggregates. 

We present the results of a test simulation, where the sticking through surface effects (S2) 
- bouncing with mass transfer (B2) transition regime in the 'pP' plot is neglected and 
replaced by S2 collisions. This alternative transition regime provides a good opportunity 
to examine the rapid growth presented in Sect. 14.4.3. H as the kernel is now simplified. 
The new kernel also provides us with the possibility to see how much the outcome of 
our simulations can be altered by changing critical areas of the parameter space. As 
the transition regime is only constrained by one experiment in a rather small area (see 
e.g. Fig. 14.51 or Fig. 11 in Chapter [3]), this part of the parameter space may require 
changing in the future. We use the same initial conditions as in the 'Mtld-4ml00' model 
described in Sect. 14.4.31 

In this case, the heaviest particle experiences increased relative velocities, as soon as it 
reaches m = 0.1 g, and the particle undergoes a rapid growth period (as in the original 
MMSN simulation, Sect. I4.4.3p . Figure 14.131 indicates the growth timescale of the 
heaviest particle (dotted line) and the minimum growth timescale possible (solid line). 
Since none of the B2 collisions can reduce the mass of the heaviest particle, the growth 
timescale reaches its possible maximum. The heaviest particle increases in mass until 
the rest of the particle population enters the B2 regimes above 0.1 g in the 'pP' plot. In 
this simulation, the maximum average mass is 27 g, whereas in the original simulation 
for the transition regime, the value was 4.18 g. 

Together with Chapter O this work is the first attempt to calculate dust growth in 
protoplanetary disks on an empirical, thus more realistic basis. However, a few additional 
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Figure 4.13: Growth timescale in the test siraulation where the S2-B2 transition 
regime is replaced by S2 cohisions only. The dotted line represents the growth timescale 
of the heaviest particles, the solid line is the minimum growth timescale. In this sce- 
nario, the growth timescale reaches its maximum possible value. 



cycles of the feedback loop between the laboratory experiments, the models of the kind 
described by Chapter O and the models described in this chapter have to be conducted 
before we can obtain a truly reliable model of dust growth in protoplanetary disks. 



4.5.2 Retention of small grains 
I I 



Dullemond & Dominik (j2005l ) showed that without a mechanism that reduces the stick- 
ing probability of particles in the upper layers of the disk or without a continuous source 
of small particles, the observed SEDs of TTauri stars should exhibit a very weak in 
frared excess. The SEDs of TTa uri stars have strong IR excess (e.g. 



Furlan et al 



(|2005l ). iKessler-Silacci et al. (|20o3)); iherefore, some kind of grain-retention mechanism 
is needed to explain these SEDs. Previous models of grain growth assumed a continu 
ous cycle of growt h and fragmentation , which provides the necessary a mount of small 
partic les (see e.g., 



Brauer et al. 



Birnstiel et al 



(|2008af ). iDullemond fc DominikI ((20051), 
(|2009l )). Our simulations, however, have shown that the mass distribution function is 
narrow. Small, monomer-sized particles are not present and fragmentation is ineffective 
in providing small particles, which could be transported to disk atmospheres. The ques- 
tion naturally arises: how can small grains be produced in our collision model? 



X^ISUSp SSDLU 'LUJOU 



One possible solution might be provided by bouncing. IWeidling et al.l (|2009l ) performed 
bouncing experiments by placing an aggregate onto an oscillating metal plate and mea- 
suring the porosity of particles due to collisions with the plate. They observed that 
approximately 10% of the projectile mass was eroded during the experiment (see Table 
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1 of their paper). This mass loss may be caused by the initial collisions, the eroded 
mass sticking to the baseplate. It is also possible that small pieces of fragments grind 
off when the aggregates bounce, which cannot be observed in the experiment. These 
ground-off particles can then diffuse out of the midplane and provide the small particles 
that we observe to enter to the upper layers of the disk. Future laboratory experiments 
are needed to quantify the level of ground-off particles created by bouncing collisions. 



The second possible explanation involves dust growth in the upper layers of the disk. 
We performed two simulations at four pressure scale-heights in both the low density 
model and the MMSN model for a = 10~^. We find that the relative velocity of two 
monomers in the Brauer model is 2 m s~^, thus monomers at these heights do not 
coagulate, only bounce. The particles in the MMSN model can form aggregates of a 
maximum of 10 /um in size. Using a higher value of a (as usually assumed in the upper 
layers of the disk), we can completely halt even this limited growth. Therefore, bouncing 
could be the key mechanism reducing the sticking probability of the particles. However, 
if substantial vertical turbulent mixing takes place, bouncing may not be able to help, 
because these monomers would then be "vacuum-cleaned" away by the larger particles 
at the interior of the disk. Additional studies of ID vertical slices of disk models are 
needed to investigate this scenario. 



4.5.3 Implications for planetesimal formation models 



It is also evident that coagulation alone is unable to produce planetesimals under the 
conditions presented in this work. Even if the turbulence parameter is assumed to be 
zero, relative velocity caused by radial drift prevents particles crossing the so called 
'meter-size barrier'. An ideal environment for particle growth is a pressure bump in the 
dead zone, where both the turbulent and radial relative velocities are reduced. This en- 



vironment is located around the snow line (jKretke &: LinI (j2007l )). iBrauer et al.l (|2008bf ) 
showed that in these pressure bumps relative velocities stayed below a fragmentation 
threshold of 10 m s~^, providing a window through which particles could overcome the 
m-size barrier, although they assumed perfect sticking (no bouncing) below the frag- 
mentation barrier. Future studies should verify whether planetesimals can be formed 
with the collision model presented in this study. 



Anot her planetesimal - formin g mechanism is the gravitational collapse of swarms of boul- 
ders (jjohansen et al. ( 2007 )). This scenario assumes that a large amount of the solid 
material is presented in dm-sized boulders {St > 0.1) at the midplane of the disk. These 
boulders then concentrate in long-lived high pressure regions in the turbulent gas and 
initial over-densities are amplified further by the streaming instability. This mechanism 
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forms 100 km sized objects on a very short timescale (some orbits). However, our sim- 
ulations produce particles with St ~ 10~^, which is due to bouncing with compaction 
(Bl) and the low (1 m s^^) fragmentation velocity of silicates. Using a 'stickier' material 
such as ices or particles wit h organic mantels may pr o duce larger particl es . Molecular 



Wada et al 



(120071) 



Wada et al 



dynam ical simulations (e.g.. iDominik TielensI (|l997l ). 

([2OO8)) have shown that icy aggregates could have fragmentation velocities of about 
10 m s~^, although these findings have yet to be confirmed by laboratory experiments. 
Similarly, it is conceivable that the enhanced sticking capabilities of ices can prevent 
bouncing, which is so omnipresent for small particles in our simulations, or shifts its 
proficiency to larger sizes. 



Cuzzi et al. (|2008l ) outlined an alternative concentration mechanism to obtain gravita- 
tionally unstable clumps of particles, which can then undergo sedimentation and form a 
'sandpile' planetesimal. In this model, turbulenc e causes dense con centrations of aero- 
dynamically size-sorted, chondrule-size particles (jCuzzi et al.l (|200ll )) or more precisely, 
particles of Stokes numbers St = Re~^^'^ ~ 10""^ in our simulations. Since growth in our 
models is typically halted at these Stokes numbers, this concentration mechanism is an 
obvious successor to coagulation - at least where it concerns the conditions adopted in 
this chapter (lAU, silicates). 



However, we emphasize that the formation of a gravitationally unstable clump does 
not imply that planetesimals will form witho ut impedime n t. An important question is 
how collisions will affect the collapse. In the ICuzzi et al.l (j2008l ) scenario, the collapse 
occurs on a sedimentation timescale and at these high densities, collisions between par- 
ticles are frequent. Likewise, in the Johansen scenario - where the collapse occurs on 
an orbital timescale and involves St ~ 0.1 particles - collisions can be rather violent. 
Collisional fragmentation or erosion may change the appearance of the collapse, because 
the small fragments are carried away by the gas. The role of collisions in these situa- 
tions is certainly an important question, and our new collision model provides a tool to 
quantitatively address this issue in future studies. 



4.5.4 Consequences for laboratory experiments 

Prom Fig. 11, in Chapter [3l one can see that only a small part of the relevant parameter 
space has been covered by experiments. Although laboratory experiments cannot be °o 
performed for every point of the parameter space, we suggest future ones based on Figs. ° 
14. 5[ 14. 7[ and 14.101 to understand dust growth in the early stages of planet formation. 



• More experiments in the 'cc' and 'cC regimes are needed as particles become 
compactified toward the end of their evolution. Thus, most of the collisions occur 
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in this regime, at velocities between 0.1 and 100 cm s , and masses between 10 
and 10 g. 

• As seen in Figs. 14.51 14.71 and 14.101 the 'hot spots', where most of the cohisions 
occur, are located in equal-sized regimes, at the left side of the fragmentation line. 
Therefore, it is important to map these areas of the parameter space in detail. 

• We define a sharp border line between the hit & stick (SI) and bouncing with 
compaction (Bl) collisions. If there is a continuous transition between SI and Bl, 
the growth of particles should not be halted by bouncing at these low particle sizes. 
Since many collisions occur in the 'pp' and 'cc' regimes, even a small probability 
of growth could increase the particle sizes. 

• As seen in Fig. 14.6b . particles in high gas density environments can have en- 
largement parameters that are much higher than 6.6 {(j) = 0.15). An interesting 
question is whether the collision types and regimes are also valid for particles with 
such low volume- filling factors, or whether these particles have different collision 
behaviors? 

• The sticking through surface effects (S2) - bouncing with mass transfer (B2) tran- 
sition regime greatly affects the outcome of the simulations (see Sect. I4.5.ip . 
However, the transition regime is only mapped at the high velocity and low mass 
regions. Therefore, it is important to constrain more tightly this part of the pa- 
rameter space. 

• The critical-mass ratio affects both particle masses and porosities. Experiments 
are needed to constrain its value. 

• The bouncing model, described in Chapter [3l has important implications for the 
evolution of dust aggregates in protoplanetary disks but is unfortunately still con- 
strained by too few experiments. Additional experiments are needed to refine the 
model, as bouncing with compaction (Bl) is the most frequent collision type in all 
of the simulations. 

4.6 Summary 

We have performed simulations of dust growth using the Monte Carlo code of ZsD08 
and a dust collision model based on laboratory experiments (Chapter [3]). We have 
performed simulations in the midplane of three disk models at low (2.4 x 10^^^ g cm~^), 
intermediate (1.4 x 10~^ g cm~'^), and high (2.7 x 10^*^ g cm^'^) gas densities at 1 AU 
distances from the central star. We have varied the turbulence parameter (a) and the 
critical-mass ratio (r^) to explore their effects on the mass and porosity distribution 
functions. Our main results are: 



118 



Chapter 4. Bouncing barrier 



• Upon using a = 10 ^, the low-density / MMSN / high-density model produces 
particles with maximum mean mass of 9.7 X g / 4.18 g / 0.23 g, respectively, 
the maximum average enlargement parameter of these particles are 7.12 / 21.9 
/ 38.0, respectively. The maximum average Stokes numbers are 2.2 x 10~^ / 
2.8 X 10-^ / 3.6 X 10-5, respectively. 

• We find that particle evolution does not follow the previously assumed growth- 
fragmentation cycles. Although catastrophic fragmentation is present for a short 
period of time in some models (typically when a = 10"^), it has a fringe effect. 
Particles in most of the simulations do not reach the fragmentation barrier because 
their growth is halted by bouncing. 

• We see long-lived, quasi-steady states in the distribution function of the aggregates 
that are caused by bouncing. The final equilibrium state is not reached within the 
drift or the viscous timescales. 

• We have performed simulations of varying turbulence strength. Wc find that the 
system is 'non-linear': the maximum mass of particles is not a decreasing function 
of the turbulence parameter and is not an increasing function of the gas density. 

• We have explored the effects of the critical mass ratio. We find that different 
critical mass ratios can aff'ect the particle evolution. Low critical mass ratios can 
produce heavier particles, while high values of Vm can halt the growth earlier. 

• The maximum Stokes number is almost independent of either the gas density or 
the strength of the turbulence. 

• The maximum mass of the aggregates is limited to 1 g because of the S2-B2 
transition regime. 

• The Stokes number 10^^ particles can be concentrated by aerodynamical size- 
sorting, thus planetesimals can form from these particles. 
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Chapter 5 

Sedimentation driven coagulation 
within the snow line 



5.1 Introduction 



Due to the vertical component of the stellar gravity, dust particles sediment towards 
the midplane of the disk. Observational evidence of the vertical sedimentation of grains 
exists for a large number of disks, although such evidence is usually indirect. 

Sub-micron grains are present at the disk surfaces as shown by scattered light images in 
the optical and near infrared (NIR) wavel engths. Such multi- w avelength scattered light 
images prov i de evi dence for grain growth (j Watson et al.l (120071 ) and references therein). 



Pinte et al 



(j2007l ) showed by reproducing multi-band images of the binary system of 
GGTau that the dust scale height for 10 micron sized particles is roughly half of that 
for micron sized particles. 



The spectral energy distribution (SED) is also affected by settling. iD'Alessio et al.l (j2006l ) 
showed that in order to explain the median SEDs of classical TTauri stars, the dust to 
gas ratio has to be reduced by a factor of 10 at the disk atmosphere compared to the 
standard value. Th ere are also indications that the settling of grains is correlated with 
the age of the disk (jSicilia-Aguilar et al.l . 120071 ). However, the connection between the 
exact shape of the 10 in i cron f eature of SEDs and sedimentation is not well understood 



(|Dullemond fc Dominikl . l2008h . 
I I 



Apai et al. (|2005l ) showed evidence for settling in disks around brown dwarf stars and 
concluded that growth, crystallization and settling of dust happens around low mass 
stars in a similar manner as around intermediate and solar mass stars suggesting that 
planet formation is a robust process. 

Sedime ntation also affects the ver tical temperature structure of the disk. The simula- 
tions of lAikawa &: Nomural (|2006l ) showed that as the dust particles sediment towards 
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the midplane, the opacity is reduced, therefore the temperature of the gas decreases. As 
the steUar radiation can now penetrate deeper in the disk, the temperature at interme- 
diate heights increases. The change in the densit y and temperature s tructure naturally 
influences the chemistry of the disk atmospheres (iBergin et al.l . 120071 ) . 



Dust sedimentation not only affects the upper layers of the disk. The gas tends to 
rotate on a sub-Keplerian velocity. The dust particles without the gas would rotate on a 
Keplerian orbit as the dust particles do not exert pressure on each other. But due to the 
presence of gas, the small dust particles move with the same (sub-Keplerian) speed as the 
gas, if the dust to gas ratio is much smaller than unity. However, if the dust to gas ratio 
approaches unity due to settling, the dust influences the motion of the gas. Therefore 
the gas at the dusty midplane layer rotates faster than the upper (not so dusty) layers 
and vertical shear is generated. This shear triggers Kelvin-Heli nholtz instability which 



Johansen et al, 



develops into turbulence. This process was firs t recognized by IWeideiischilling 
and it is still an actively researched area (e.g. by 



( 200' 



Chiand(200 



(Il980h 



In the framework of this thesis, a collaboration started between the lab-community 
and the modelers to better constrain the dust evolution in protoplanetary disks using 
a realistic collision model tha t is based on the laboratory experiments. In Chapter [3] 
(based_on 



on 



Giitt 



Zsom et al, 



er et al 



(|2010l )). we introduced this collision model. In Chapter H] (based 
(|20ld )) , we used this collis i on in odel for the first time in the Monte 
Carlo (MC) method of lzsom fc Dullemondl ([2003) (henceforth ZsD08). The models of 
Chapter H] were local box models meaning that the dust evolution was only followed at 
one location of the disk. These models showed that bouncing plays an important role 
in dust evolution. 

We further develop these models to simulate a ID vertical column in the disk thus 
investigating sedimentation driven coagulation. We want to better understand the pro- 
cess of sedimentation and the role of particular physical phenomena like porosity of the 
aggregates, collision models and turbulence. 



Previous work by 



Dullemond &: Dominik 



(|2005f ) (henceforth DD05) showed that without 
a mechanism that reduces the sticking probability of particles in the upper layers of the 
disk, or without a continuous source of small particles, the observed spectral energy 
distributions (SED) of TTauri stars should exhibit a very weak I R excess. In c o ntrast , 
the observed SEDs of TTa uri stars have strong IR excess (e.g. iFurlan et al.l (j2005l ): 



Kessler-Silacci et al 



20061 )) therefore some a grain-retention mechanism is needed to 



explain the SEDs. 

Previous models of grain evolution assumed a continuous cycle of growth and fragmenta - 
tion, which pro v ides t he necessary amount of small particles (e.g. 



Brauer et al. 



Birnstiel et al. 



(2008a 



(I2OO9I )). However, Chapter [3] and [H showed that particle evolution is 
halted by bouncing and no cycle of growth and fragmentation is present. We simu- 
late dust evolution driven by Brownian motion, turbulence, and sedimentation in a ID 
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vertical column of the inner disk. We investigate the time evolution of sedimentation 
driven coagulation, and search for ways that can keep a sufficient amount of the small 
dust particles levitated at several pressure scale heights to explain the observed SEDs 
of young stars. 

The chapter is organized as follows. In Sec. 15.21 we describe the numerical method used 
to follow the particle motion and coagulation. We validate the code and increase the 
complexity of the model step-by-step in Sec. 15.3.11 We show the results in Sec. 15.31 
finally we discuss those results in Sec. 15.41 and provide a summary in Sec. 15.51 

5.2 Numerical method 
5.2.1 Basic considerations 

The local box approach in Chapter [3] is based on two assumptions. 1.) The particles 
are homogeneously mixed inside the box. 2.) Particles do not enter or leave the box, 
i.e. it is closed. Due to these two assumptions, it was not necessary to follow the exact 
location of the particles. 

In the models considered here, however, we place such boxes (or grid cells) on top of 
each other to simulate a ID column in the disk and follow how particles settle towards 
the midplane. Inevitably, particles move from box to box during this process. Therefore 
the assumption that particles cannot enter or leave the boxes has to be relaxed. The 
first assumption of the method in Chapter [3] is kept, we still assume that the particles 
inside a given box are homogeneously distributed when we consider coagulation (for 
sedimentation, the individual positions of the particles are used). The second assumption 
is modified in the following way. 2.) The simulated column is closed, e.g. particles inside 
the column can move freely vertically. However neither do new particles enter from the 
"outside" , nor do particles from inside the column leave. As particles move through the 
boxes, it is necessary to follow the position of the particles (see Sec. I5.2.3P as we must 
find out in which box a particle is located. 

The motion of particles imposes a limit on the time-step of the simulation. We do not 
want the particles to move more than one box in a time-step. A sedimenting particle 
should have the possibility to interact with all other particles along its way, it should 
not skip over boxes thus avoiding the particles in it. Therefore, we use an adaptive time- 
stepping method. The maximum of all particle velocities is obtained (t'max); and since °° 
we know the height of the boxes (/ibox); the maximal (safe) time-step can be determined - 



as 




At = C 



^^max 



where C is the Courant number which we typically set to be 0.1. 
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The code schematically performs the following steps: 

1. First the velocities of the particles are calculated. 

2. A safe time-step is determined to avoid particle 'jumps'. 

3. The position of the particles is updated using their velocities, their previous posi- 
tions, and the time-step. 

4. We determine the box in which each particle resides. 

5. We call the coagulation subroutine described in Chapter H] to calculate the evolu- 
tion of the particles separately in each box for the given time-step. 



5.2.2 Initial conditions 



We assume that the gas density profile is constant during the simulation. This assump- 
tion is valid if the simulated time is less or comparable to the viscous timescale of the 
gas. The viscous timescale can be calculated as 

tvis = r^UT, (5.2) 

where r is the distance from the central star, ut is the turbulent viscosity. A typical 
val ue for tyjs at 1 AU is lOf - lO'^ yrs. We assume that turbulence is parameterized by 
the IShakura fc SunvaevI (jl973l ) a parameter 

i^T = aCsHg, (5.3) 

where Cg is the isothermal sound speed, and Hg is the pressure scale height of the gas 
disk. The turbulence parameter a reflects the strength of the turbulence in the disk. 
Typical values range between a = 10~^ and 10~^, where the former corresponds to the 
turbulent strength in dead zones, the latter describes turbulence in disk atmospheres. 

The vertical structure of the disk is determined by the equilibrium between the vertical 
component of the gravitational force and the acceleration due to the vertical pressure 
gradient in the gas. If the disk mass (Mjisk) is much smaller than the mass of the 
star (M*), and the vertical thickness of the disk (Hg) is a small fraction of the radial 
distance (both conditions are safely met for the disk parameters described below), then 
the vertical density can be approximated as 

Pg{r, z) = -^l- eM-^V-^H^g), (5.4) 

where S(r) is the gas surface density at distance r, and z is the height above the 
midplane. In this chapter we choose M* = 0.5Mq, r = 1 AU, S(l AU) = 100 g/cm^ 
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similarly to DD05. The pressure scale height can be calculated as 

Hg = Cs/n, (5.5) 
where ri is the orbital frequency at 1 AU. The isothermal sound speed is 

Cs = J-^, (5-6) 

where kB is the Boltzmann constant, /i is the molecular weight, which is 2.3 for molecular 
gas, nip is the mass of the proton, and T is the temperature of the gas, which is 200 
K for the stellar and disk parameters considered above (see DD05). We assume the 
temperature to be constant as a function of height. This is a reasonable assumption if 
the temperature of the gas is solely determined by the stellar irradiation. 

We simulate 4 pressure scale heights (0.16 AU above the midplane), use 40 evenly spaced 
boxes and 10^ particles in the simulations unless otherwise stated. The number of boxes 
and the number of particles are chosen by taking into account two points. 1.) The upper 
box is initially not empty, it contains at least a few particles. 2.) Simulations, which 
are performed with the exact same initial set-up for the gas and using the same collision 
model, but using different initial positions for the dust (e.g. using a different seed for 
the random number generator), do not differ qualitatively. Differences are expected 
due to the intrinsic stochasticity of the Monte Carlo implementation. Every particle 
represents the same portion of the total dust mass, therefore more particles are present 
in the lowest box and only a few in the upper box. The initial particle positions are 
determined randomly, therefore the initial gas to dust ratio is noisy. The noise is lower 
closer to the midplane, and gradually increases with height. However, the mean initial 
gas to dust ratio is constant, 1:100. 



5.2.3 Position update 

The position of the particles are determined by vertical settling and turbulent diffusion. 
In principle, Brownian motion also contributes to the change of particle positions, but 
its effect is negligible compared to the other two effects. 



The equation g overning the diffusion , and settling of the dust in a n on-homogenous gas 



density field is (jPubrulle et al 



1995 



Fromang &: Papaloizoul . l2006l ) 



dpd ^ d 
dt dz 

or in a more practical form: 
d^DdPd 



DdPr 



d ( Pd 



dz 



Pg 



d 

+ —{^ihsPdz) 



dpd 
dt 



dz^ 



dz \ Pg dz 
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where pd is the dust density, 0^ is the diffusion coefficient of the dust (for the calculation 
of Dd, see the next paragraph) and ts is the stopping time of the particle. The stopping 
time is the timescale a particle needs to react to the changes of the surrounding gas. We 
define the dimensionless Knudsen number being 

Kn = (5.9) 

where a is the size of the aggregate, and Amfp is the mean free path of the gas. A particle 
is in the Epstein regim e if Kn < 1 (to be more precise, if a < |Amfp), where the stopping 



time is (jEpsteinl (|l924l )): 

where m and A are the mass and the aerodynamical cross-section of the particle, and Uth 
is the thermal velocity. If the Knudsen number is greater than 1 (at high gas densities 
where the mean free path is low or in the case of large particles), the first Stokes regime 
applies and the stopping time becomes 

3m 4 a , x 

ts = tst = -. jx-- . (5.11) 

4vthPgA 9 Amfp 



The first term on the right hand side of Eq. 15.81 is the well-known diffusion term. Using 
only this term, particles with tg = (tracers) would be homogeneously distributed as 
a function of height over several diffusion timescales. The first and the second term 
together on the right hand side ensures that the tracer particles will be distributed 
according to the background gas density field. The third term describes the settling of 
the particles. Equation l5.8l is valid if the motion of the dust does not influence the motion 
of the gas (the back-reaction from the dust to the gas is negligible). This condition is 
met if the dust to gas ratio is ^ 1. 

We note that we do not solve for this dust density field directly. We follow the motion 
of dust particles, each of which represents a portion of the total dust mass inside the 
column, thus we derive the corresponding velocities (or fluxes) for the first, second, and 
third terms of Eq. 15.81 to calculate the position update of the particles. 

We calculate D^, the diffusion coefficient of the dust, an d define the diffusion velocity , 
V£)i. The diffusion coefficient of the gas can be defined as (iDullemond &: Dominikl . l2004l ) 



Dq = h'T = aCsHq 



(5.12) 



Based on 
as 



Youdin &: Lithwick 



(j2007l ) , the diffusion coefficient of the dust can be calculated 
Dd = Dg/{l + St'), (5.13) 



where St is the Stokes number 



St = ts^. 
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The average displacement of a particle in ID during the time-step of At then is 



L = y/2DdAt. (5.15) 

The real displacement of the particle (Az) is drawn from a Gaussian distribution which 
has zero mean and a half width of L. The "diffusion velocity" can then be calculated as 

VDi = Az/At. (5.16) 

This velocity component tries to smear out dust concentrations. It is important to note 
that the real, physical velocity of the particle during turbulent diffusion changes ran- 
domly every time the aggregate interacts with a turbulent eddy. The diffusion velocity 
defined above is a numerical construct to calculate the time-averaged velocity of the 
particle during a time-interval At. 

The second term on the right side of Eq. 15.81 results in a systematic velocity term which 
pushes particles towards the density maxima of the gas. This velocity can be determined 
by 



VD2 = Dd—^- (5.17) 

Using these two velocity components {vdi and vd2)^ the particles with non-zero stopping 
times will be distributed according to the gas density profile in a timescale longer than 
the diffusion timescale. 

The fact that particles with > settle towards the midplane and have a scale height 
less than Hg is the result of the third term of Eq. 15. 8^ the settling velocity. The settling 
velocity of a particle can be determined by 

7;set = -z^hs. (5.18) 

The new position of the particles can then be determined by using these three velocity 
terms: 

Z = Zo\d + {VDI + VD2 + ^^set)At. (5.19) 



5.2.4 Coagulation 

The collision model used in this work is similar to the one used in Chapter HI There are 

X^jSUSp SSDLU 'LUJOU 

however two differences. The first difference is the additional source of relative velocity % 'a b L o 
due to differential vertical settling (see Eq. I5.18p : 5 2 2 5 5 



Avs = \Vsetl - Vset2\- (5.20) 



One could use the individual height of the particles (as calculated in Sec. I5.2.3P to obtain 
Avs, but that would violate the first assumption of the coagulation model (namely that 
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the particles are uniformly mixed inside the box). When particles collide, they have 
a zero distance, therefore their height must be identical. If the (non-representative) 
particles are uniformly mixed, two colliding particles with identical stopping times do 
not have differential settling velocity. If, however, one uses the individual heights of the 
representative particles, it is implied that the particles are not uniformly mixed in the 
box. Their height can be different which would result in a non-zero settling relative 
velocity if tgi = ts2- Therefore, we must use an averaged height, the height of the box, 
to reliably calculate Avs- 

The second difference is in the calculation of the aerodynamical cross section which is 
used to calcu late the stopping ti me. In Chapter H] we used geometrical cross section of 
the particles (jormel et al. . 2007 ) 



A 



(5.21) 



where Tc is the compact radius of the aggregate (assuming that the mass of the particle is 
contaii ied in a compact sp here of radius Tc), and ^ is the enlargement parameter defined 

dimension 



m e.g. 



Ormeletal 



(120071 ). This formula works well for par ticles with fracta . 



above 2. However, we will also use the porosity model of lQkuzumi et al.l (|2009l ). and 
their model produces aggregates with fractal dimension below 2. If one calculates the 
stopping time of such an aggregate in the Epstein regime (Eq. IS.lOp using the formula in 
Eq. 15. 2H one gets that the stopping time is less than the stopping time of a monomer. 
This is clearly unphysical. The reason for this low stopping time (large area) is that 
Eq. 15.211 does not take into account the empty space between the 'fractal branches'. 
To avoid such unphysical results for aggregates wit h low fractal dimensio ns, we use the 
aerodynamical cross section as defined in Eq. 47 in 



Qkuzumi et al 



(120091). 



5.3 Results 



We perform 11 simulations, in which we gradually use more realistic collision models, 
investigate the effects of different porosity models and turbulence. The IDs and param- 
eters of these simulations are shown in Tab. 15.11 First we compare our model against 
the r esults of DD05 (model DP in Tab. 15.10. Then we use the porositv model of Ormel 



et al. 120071 ) and 



Qkuzumi et al 



(|2009l ) to investigate the effects of porosity (models 
DDa and DDb, respectively). So far we assume that the aggregates stick together at all 
relative velocities and the turbulence parameter a is zero. In the next step we construct 
a more realistic collision model with sticking, bouncing and fragmentation (models SBl, 
SB2). We call this collision model the "simplified Braunschweig model" because it uses 
only three collision types out of 9, which is described in Chapter [3] (the complete Braun- 
schweig model). In the next step we turn on turbulence (models SB3-6) to examine 
the effects of turbulent stirring. Finally we use the complete Braunschweig model with 
turbulence (models FBI, FB2). 
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Table 5.1: Overview and results of all the sedimentation simulations. 



ID 


Coll. model 


Por. model 


a 




^rain 


^rain 




Hp 












[yrs] 


[g] 




[Hg] 


(1) 


(2) 


(3) 


(4) 




(5) 


(6) 


(7) 


(8) 


DD 


hit&stick 


— 







500 


Ti T 

10-^ 


1 





DDa 


hit&stick 


Ormel 







600 


10" 


10 





DDb 


hit&stick 


Okuzumi 







900 


10^ 


106 





SBl 


simpl. Br. 


Ornicl 







3000 


10-5 


50 





SB2 


simpl. Br. 


Okuzumi 







6000 


10-2 


10^ 





SB3 


simpl. Br. 


Ormcl 


10- 


-6 


2000 


10-5 


30 


0.1 


SBl 


simpl. Br. 


Okuzumi 


10- 


-6 


4000 


10-2 


10^ 


0.1 


SB5 


simpl. Br. 


Ormel 


10- 


-4 


500 


10-6 


10 


0.5 


SB6 


simpl. Br. 


Okuzumi 


10- 


-4 


700 


10-4 


10^ 


0.5 


FBI 


compl. Br. 


Ormel 


10- 


-4 


500 


10-6 


10 


0.25 


FB2 


compl. Br. 


Okuzumi 


10- 


-4 


700 


10-4 


10^ 


0.25 



Col. 1 is the ID of the simulations, col. 2 describes the used collision model (hit&stick, 
simplified Braunschweig model, or complete Braun sch weig model), col. 3 indicates the 



used porosity model (based on Ormel et al. ( 2007 ) or Okuzumi et al. ( 20091 )). col. 4 



describes the value of the turbulence parameter a, col. 5 is the time the rain-out particle 
reaches the midplane of the disk in years, col. 6 and 7 are the mass and the enlargement 
parameter of the rain-out particles, respectively, and col. 8 is the scale height of the dust 
expressed in the scale height of the gas. Note that for the FBI and FB2 simulations, 
the 

i^vain values are not the final masses of the particles, but the masses the rain-out 
particles have when they first reach the midplane. The final masses of the particles are 
10-^ g for both the FBI and FB2 simulations. 



5.3.1 Test: comparison with the DD05 model 



DD05 performed a simulation (S2 in their paper, DD in this Chapter), where the disk 
model is the same as the one described in Sec. 15.2.21 the particles are compact, upon 
collision particles stick together at all collision energies, and the only source of relative 
velocity is Brownian motion and differential settling. They found that the 'rain-out' 
particles reach the midplane in 500 yrs having attained sizes of a millimeter (10~^ g in 
mass) . 

We performed the exact same simulation to validate our code. We find that the 'rain- 
out' particles in our simulation reach the midplane also at 500 yrs and have masses of 
10~^ g, therefore we can conclude that our code works properly. We illustrate the mass 
evolution of particles as a function of their height at six different snapshots {t = 1 yr, 
10 yrs, 100 yrs, 316 yrs, 10^ yrs, 10^ yrs) in Fig. I5.1[ 

Brownian motion is essential in our simulations because growth by Brownian motion 
kicks in the sedimentation driven coagulation. The reason is that we have initially a 
mono- disperse particle size-distribution (meaning that all monomers have the same size 
and mass), therefore the aerodynamical properties of the monomers are identical, thus 
there is no relative velocity due to settling between the monomers at a given height. 
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Figure 5.1: The mass distribution at i = 1, 10, 100, 316, 10^ and 10" yrs for the DD05 
model. The x axis is the mass of the aggregates in grams, the y axis is the height above 
the midplane expressed in units of the pressure scale-height. The contours represent 
the normalized mass density of the dust. The sub-figures illustrate the dust to gas ratio 
(x axis) as a function of height above the midplane (y axis) 



If growth due to Brownian motion was not initiated (e.g. growth by Brownian motion 
did not introduce aggregates with different aerodynamical properties than that of the 
monomers), the monomers would simply sediment to the midplane without any growth. 
Although DD05 included Brownian motion, this effect would not have been present as 
that simulation started with a (narrow, but not infinitely narrow) size distribution. 

As shown in Fig. 15. 1^ growth by Brownian motion is faster at the midplane due to the 
higher gas and dust densities (at t = 1 and 10 yrs). Once particles at the upper layer also 
start growing by Brownian motion, sedimentation driven coagulation starts and particles 
at the upper layers grow much faster than the aggregates at the midplane (at t = 100 
yrs). The heaviest particles sweep up the smaller particles while they sediment and 
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further increase their setthng velocity, resulting in a rain-out at t = 500 yrs. Once the 
first rain-out particles reach the midplane, they could only grow by Brownian motion, 
because at the midplane, the settling velocity of any particle is zero (see Eq. I5.18p . But 
the relative velocity due to Brownian motion for such heavy aggregates is low, therefore 
the particles that have reached the midplane, do not increase in mass significantly during 
the simulation. 



5.3.2 The effects of porosity 

In the previous section we used compact particles. However fluffy particles couple better 
to the gas. Therefore we perform simulations with two different porosity models to 
investigate the effects of porosity. We include growth by Brownian motion and settling 
only, and assume that particles stick together at all collisional energies. 



The Ormel model - DDa First we use the porosity model of lOrmel et al.l (j2007l ). 

This porosity model is based on both theoretical and experimental investigations of 
the microphysics of dust aggregates. It is essentially incorporates PCA-like collisions 
(particle-cluster aggregation - collisions between particles and clusters) and CCA-like 
collisions (cluster-cluster aggregation - collisions between clusters of similar size) and 
interpolates between these two. In this porosity model, compaction can occur if the 
collision energy is Ecoii > 5£Jroih where E'roU is the rolling energy, which is the energy 
needed to roll two monomers by 90 d eg. However, we do no t include this regime to stay 



consistent with the porosity model of lOkuzumi et al.l (120091 ) . as their model also do not 
treat compaction. 

The evolution of the mass can be seen in Fig. 15.21 the evolution of the enlargement 
parameter is illustrated in Fig. 15.31 The evolution of these porous dust particles shows 
a similar behavior to the evolution of compact particles. The rain-out happens somewhat 
later at t = 600 yrs, and the particles that reach the midplane have masses of 1 g, two 
orders of magnitud e higher than in the previous model, which agrees with the results of 



Ormel et al 



(120071). 



The porosity evolution (Fig. 15. 3p shows some interesting features. The enlargement 
parameter naturally increases during Brownian motion (at t = 1, 10 yrs) and also 
during the initial phases of settling (at t = 100 yrs). However, as the particles approach 
the midplane, they become more compact (at t = 10^ yrs). The particle population that 
the rain-out particles can sweep up is getting smaller as they approach the midplane 
(see Fig. 15.21 at t = 10^ yrs). Therefore the rain-out particles collide with ever smaller 
particles and these small particles can 'fill up' the holes of the rain-out particles, the 
collisions are more PCA-like and the enlargement parameter decreases. 
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Figure 5.2: The mass distribution for the DD05 model using the Ormel porosity 
model. The axes and contours are the same as in Fig. 15.11 



The Qkuz umi model - DDb The second porosity model we use is the one con- 
structed by 



Qkuzumi et al 



(|2009l ). This collision model defines a third type of ag- 
gregation next to the PCA and CCA collisions, that is QCCA (quasi cluster-cluster 
aggregation - collisions between clusters with a predefined mass ratio). This model is 
based on numerical models of geometrical sticking (e.g. no restructuring of monomers 
happens during a collision). 

The mass and the porosity evolution are shown in Figs. 15.41 and 15.51 The most striking 
property of this simulation is the maximum mass and porosity of the particles. We 
end up with particles of 10^*^ g in mass having an enlargement parameter of almost 
10^ (the compact radius of such an aggregate is some meters, however, the enlarged 
radius is several kilometers). The stopping time in the Epstein regime (Eq. I5.10p 
is proportional to m/A. As the Okuzumi-model produces aggregates with a fractal 
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enlargement parameter enlargement parameter 

Figure 5.3: Enlargement parameter distribution for the DD05 model using the Orniel 
porosity model. The x axis here represents the enlargement parameter of the aggregates, 
and the contours show the normalized mass-weighted enlargement parameter of the 

particles. 



dimension of ~ 2 (the mass scales with a^, where o is the particle radius), the stopping 
time only slightly increases with mass. Therefore particles settle slowly and produce 
extremely fluffy structures. At one point, however, the particle size becomes larger than 
the mean free path of the gas, and the aggregate enters the Stokes regime (Eq. I5.11|) . 
As the stopping time is now proportional to ma/A, the stopping time more strongly 
increases with mass. This transition from Epstein to Stokes regimes happens at t=900 
yrs for the particles located at 1.7 Hg above the midplane. Once the transition happens 
for a given particle, it settles to the midplane in a matter of years due to the heavy mass 
of the aggregate. Therefore a = Amfp is a natural upper size limit to expect rain-out in 
this model. Such huge particles are in reality probably very fragile and would break up 
due to the smallest perturbations in the gas (e.g. turbulent eddies, or the "fall" to the 
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Figure 5.4: Mass distribution for the DD05 model using the Okuzumi porosity model. 
The axes and contours are the same as in Fig. 15.11 



midplane after entering the Stokes regime). 

If one uses only the Epstein stopping time (although this is not physical), the particles 
do not sediment so rapidly to the midplane, their fractal growth goes on unhindered. 

Although porosity can somewhat delay sedimentation (the particles in the DDa and DDb 
simulations reach the miplane later than in the DD simulation), this delay is limited. A 
monomer has the lowest available stopping time, as it has the lowest settling velocity. 
In the disk model used in this chapter, a monomer needs 1.5 x 10^ yrs to settle from 4 
Hg to 2 Hg. No matter what porosity model we use, the sedimentation timescale down 
to 2 Hg cannot be longer than this. 

We must emphasize that any collision model containing exclusively sticking is only valid, 
if no significant restructuring happens during the collisions (e.g. the collision energy is 
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10° 10^ lO'* lO*" 10*^ 10° 10^ 10^ 10*^ 10*^ 



enlargement parameter enlargement parameter 

Figure 5.5: Enlargement parameter distribution for the DD05 model using the 
Okuzumi porosity model. The axes and contours are the same as in Fig. 15.31 



less than SE'rolh the rolling energy, which is the energy needed to roll two monomers by 
90 deg). This condition is clearly not met at all times in our simulations, e.g. the rain- 
out particles can have collision velocities with the swept up particles as high as several 
10 m/s in these simulations. Such collisions would result in catastrophic fragmentation. 
Therefore the results presented in this section should be considered as toy models. 



5.3.3 A simplified Braunschweig model 



X^ISUSp SSDLU 'LUJOU 



In this section we construct a simplified version of the collision model described in 
Chapter [3l We assume sticking, if the collision energy is smaller than 5ii^roii- Bouncing 
with compaction is used if the collision energy is greater than S-Eroih but the relative 
velocity of the two aggregates is less than 1 m/s. Fragmentation occurs if the relative 
velocity of two aggregates is greater than 1 m/s. The recipe for mass and porosity 
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Figure 5.6: Mass distribution for the simplified Braunschweig model using the Ormel 
porosity. The axes and contours are similar to Fig. 15. f I Notice the x axis ranges from 

10^^'^ g until 10 g in this figure. 



evolution for bouncing and fragmentation is taken from Chapter O (our hit & stick 
(SI), bouncing with compaction (Bl), and fragmentation (Fl) collision types). We still 
assume that the particles grow by Brownian motion and settling only (the effects of 
turbulence is discussed in the next section), but we use both porosity models discussed 
in the previous section. 



The Ormel model - SBl The evolution of the mass in case of the Ormel porosity 
model is shown in Fig. 15.61 The mass distributions at t = 1, 10, 100 yrs are identical 
to Fig. I5.2i At t = 316 and 1000 yrs we see the effects of bouncing at the intermediate 
energies. The rain-out particles cannot increase their mass, if they suffer bouncing 
collisions. Therefore the rain-out particles only reach masses of 10~^ g when they arrive 
at the midplane. As the rain-out particles are smaller, they settle slower, therefore 
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Figure 5.7: Mass distribution for the simplified Braunschweig model using the Ormel 
porosity with a = 10""*. The axes and contours are the same as in Fig. 



they reach the midplane only at t = 3000 yrs. The enlargement parameter is affected 
by bouncing, but the typical enlargement parameter is between 10 and 100 as in the 
DD05a model in Fig. [Ol 



The Okuzumi model - SB2 The Okuzumi model produces fiuffier particles. As a 
result the particles are heavier (10~^ g) but arrive to the midplane later (at t = 6000 
yrs) than in the previous simulation. 
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5.3.4 The effects of turbulence 
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So far all particles sooner or later ended up at the midplane because there was no effect . mb,BM 
that could counteract settling. In this section we examine the effects of a non-zero 
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Figure 5.8: Mass distribution for the complete Braunschweig model using the Ormel 
porosity with a = 10~^. The axes and contours are the same as in Fig. 15.61 The 
snapshots of this simulation is presented at the lower right corners of the thesis as a 

flip-cartoon. 

turbulence parameter, which can stir particles back up. 

A small turbulence parameter (a = 10^^) does not affect significantly the masses of the 
rain-out particles compared to models SB 1-2 of the previous section (see models SB3-4 in 
Tab. 15. ip . As particles do not only settle but also diffuse downward (and upward) due to 
turbulence, the time the rain-out particles reach the midplane is somewhat shorter than 
for models SBl and 2. In all previous simulations an infinitely dense dust layer formed 
at the midplane of the disk. However, even this low level of turbulence can prevent the 
formation of this layer and introduce a non-zero (although small) dust scale-height. 

The influence of turbulence is more pronounced if a = 10~^. The mass evolution of 
the aggregates using the Ormel porosity model is shown in Fig. 15.71 The first rain-out 
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particles reach the midplane already at t = 500 yrs due to downward diffusion, although 
these particles have lower masses than in model SBl (10~^ g - therefore, in the absence 
of turbulence, these particles would reach the midplane later than the particles in SBl). 
The dust distribution reaches a steady state at t = 10*^ yrs. 

We see that the particle mass is constant as a function of height at t = 10^ yrs. As 
turbulence effectively mixes the particles, and as bouncing prevents further growth or 
fragmentation (the dust growth is halted), both the masses and porosities of the aggre- 
gates are similar at all heights. This has an important consequence for observations. If 
the t urbulence param eter is constant as a function of height (which might not be true 



see 



Gammid (j 19961 ) and Sec. I5.4.2p . we expect that the particles observed at the disk 
atmosphere have the same properties as the ones located at the disk midplane. If how- 
ever a is some function of the height (e.g. high at the disk atmosphere and low at the 
dead zone), we might not be able to constrain the particle properties at the midplane of 
the disk, unless the disk is optically thin at the given wavelength. 

We also see from these simulations that a higher turbulence value reduces the mass of 
particles and increases the dust scale height. If turbulence is strong enough (a ~ 1), the 
dust scale height can be similar to the gas scale height and the disk atmosphere remains 
dusty at all times. However, such high turbulence value prevents any significant dust 
growth, which is not a fertile environment for planet formation. 



5.3.5 The complete Braunschweig colhsion model 

In this Section we use the complete Braunschweig model (see Chapter [3] for details), the 
value of the turbulence parameter is a = 10""^ and calculations are performed with both 
the Ormel porosity model (FBI) and the Okuzumi porosity model (FB2). 

In the simplified Braunschweig collision model, the growth is halted by bouncing immedi- 
ately if the particles enter the bouncing regime. However, in the complete Braunschweig 
model, there are possible ways for growth beyond the hit&stick border line (that is where 
-E'coii > S-Ej-oii). The most important area is in the 'pP' regime where a small porous 
projectile collides with a heavy porous target (see Fig. I3.1ip . Due to these "green" 
areas at intermediate collision energies, particles in the FBI and FB2 simulations grow 
to higher masses than in the SB5 and SB6 simulations. As a consequence, the scale 
height of the dust is lower in these simulations, as heavier particles are more difficult to 
stir up by turbulence. We illustrate the mass distribution at t = 1, 10, 100, 316, 10^, 
10"^ yrs in Fig. ESI 

The particle evolution has two phases in these simulations. The ffist 1000 yrs are 
identical for the SB5/SB6 and FB1/FB2 simulations, respectively (see also the ffist five 
snapshots of Figs. 15.71 and f5.8|) . In this phase, particles start sedimenting, and the rain- 
out particles reach the midplane. The dust evolution in the SB5 and SB6 models halts at 
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Figure 5.9: The collision history for the FBI simulation. The eight regimes of the 
complete Braunschweig model are shown with the corresponding border lines of the 
nine different collision regimes (white solid lines) . The x axis is the relative velocity of 
the particles in cm/s, the y axis is the mass of the projectile in gram units. The grey 
boxes indicate the areas that are covered with laboratory experiements (see Chapter [3] 
for more details). The colors indicate how many collisions happened at the given part 
of the parameter space during the simulation. The red and yellow areas are 'hot spots', 
where most of the collisions take place. 
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Figure 5.10: The colhsion frequency of the nine cohision types. The x axis is the time 
in years, the y axis indicates the cohision types. The colors of the stripes indicate the 
cohision frequency (e.g. the number of cohisions per year). The cohision frequency is 
shown at the midplane (a.), at 1 pressure scale height above the midplane (b.), and at 

2 Hg (c). 



this point as only bouncing collisions happen. However, during the second phase of the 
FBI and FB2 simulations, a short run-away growth appears as the relative velocities are 
"boosted" at this point (see Sec. 14.4.3.11 for a detailed explanation). Due to the rapid 
growth, particles reach 10~^ g in mass for the FBI and FB2 simulations. 

Figure 15.91 illustrates the collision history of the FBI simulations. If we compare this 
figure to Figs. 14.5^ 14.71 and 14.101 we see that the features are more smeared out in Fig. 
15 .91 than in the other figures. As we simulate here several boxes at different heights above 
the midplane, the physical conditions (e.g. gas density and sedimentation velocity) at 
the midplane and at the upper scale heights of the disk are different, which is responsible 
for the smeared out features of Fig. 
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It is more interesting to investigate the colHsion frequency of the nine cohision types as 
a function of time (see Fig. IS.lOp . If one compares this figure with Figs. 14.4b . 14.6b . and 
14.9b . we immediately see that the diversity of occurring colhsion types is much greater in 
the FBI (and also in the FB2) simulation, although the strength of the turbulence is the 
same in all cases (a = 10~^). This can be explained by the presence of sedimentation. 
The particle population in a given box is not fixed as in Chapter HI During the rain-out 
process, heavy particles coming from 1-2 Hg 'hit' the particle population at the midplane. 
Due to this process, the relative velocity between the small, midplane particle population 
and the generally larger rain-out population is increased, thus collision types can occur 
that require larger collision energies. 



5.4 Discussion and future work 



5.4.1 Porosity models 



We examined the effects of two porosity models on settling in this Chapter. Here we 
critically discuss the strong and weak points of these models. 



The porosity model of lOrmel et al.l (|2007l ) only treats PCA and CCA collisions (particle- 
cluster and cluster-cluster aggregat ion, respectively) and "constructed" semi-analytical 
recipes for intermediate size ratios. lOkuzumi et al.l (j2009l ) improved on this by modeling 



the quasi-CCA collisions (QCCA), where two clusters with a given mass ratio can collide. 
QCCA growth generally produces fluffier particles. 



However, the porosity model of lOkuzumi et al.l (j2009l ) assumes collisions of pure geo- 
metrical sticking (no restructuring) . This assumption is correct if the collision energy is 
always much lower than the rolling energy, e.g. during particle growth by only Brownian 
motion. In that case the collision velocity of the aggregates is decreasing as the particles 
grow, which ensures that the collision energy is always much smaller than the rolling 
energy. However, if particles grow also by differential settling, turbulence, or differential 
radial drift, the assumption might not be correct. If for example a collision happens with 
an energy just slightly smaller than 5£'roii, the restructuring for a single collision can be 
negligible. However, if an aggregate experiences many such collisions, the cumulative 
effect of the small restructuring might not be negligible anymore, as the restructuring 
does not behave as a random error, but it is additive in nature. 



The porosity model of lOrmel et al.l (|2007l ) is more robust because restructuring at 
energies ~ -EroU can be included (but we do not include it, see Sec. I5.3.2p . however it 
most probably underestimates the enlargement parameter by not taking into account 
the QCCA collisions. 

We propose the development of a porosity model that combines the advantages of both 
models, e.g. it incorporates QCCA collisions, but does not assume geometrical sticking. 
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5.4.2 Where can a be constant as a function of height? 
I I 



Gammie (| 19961 ) proposed the concept of layered accretion disks. If the ionization fraction 
of the gas is not sufficient to support magneto-rotational instabihtv fMRI - Balbus &: 



Hawley (|l99ll )). the turbulence parameter drops and a dead-zone forms at the midplane 
of the disk. The extent of the dead-zone is uncertain, as the ionization processes of the 
gas are not well-constr ained. For typical TTauri disks it can extend between 0.1-4 AU 



(D'Alessio et al 



19981 ). Inside 0.1 AU, the thermal radiation from the star can keep the 



dust sufficiently ionized for MRI, and outside 4 AU, the gas surface density is typically 
below 100 g/cm^, therefore cosmic rays can penetrate the disk and keep it MRI active 
at all heights. 

If a disk can be observationally resolved down to 4 AU or inside 0.1 AU and one can be 
sure that no dead-zone is present at the resolved location of the disk, the dust properties 
at the whole column can be constrained even if the disk is not optically thin at the given 
wavelength (see Sec. I5.3.4p . 

It would be however interesting to investigate how dust evolves in a layered disk model. 
Small dust p a rticles can very efficiently sweep up charges in the gas. As shown by 



Turner et al. 



(|2010l ). the dead-zone can extend to 2 Hg for 1 micron sized particles, 
but it shrinks below 0.5 Hg for aggregates that are 100 micron in size. In a simulation 
like the one presented here, this would mean that as the particles grow, the dead-zone 
shrinks. When the dead-zone disappears, the whole disk becomes MRI active and the 
particles settled to the midplane might be fragmented and stirred back up. This could 
lead to an oscillatory process. 



5.4.3 Can bouncing keep the disk atmospheres dusty? 



Another way to keep disk atmospheres dusty might be possible via bouncing. In the 
previous models of dust evolution, the c ycle between g r owth a n d fragmentation p rovidi ed 
a source for small particles (see e.g. 



Brauer et al 



(|2008a|); 



Birnstiel et al 



(120091)) 



However, the dominant collision type in our simulations is bouncing. 
I I 



Weidling et al. (|2009l ) observed that the mass of the particles in their multiple bouncing 
experiments were reduced by the end of the experiment and it was unclear why this 
happened. We propose that small pieces can grind off from the aggregates while they 
bounce. These grind off bits and pieces can provide a continuous source of small particles. 

It is possible to theoretically constrain the critical amount of grind off particles that is 
necessary to keep the disk atmospheres dusty using the models presented in this Chapter. 
In more detailed laboratory experiments, it should be possible to measure whether that 
mass- loss is really present. 
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5.5 Summary 



We performed simulations in a ID vertical column of a protoplanetary disk to better un- 
derstand the process of sedimentation. We simultaneously solved for the particle motion 
and growth inside this column. The complexity of the models was gradually increased 
to examine the effects of different processes. The first simulation used a collision model 
that only contained sticking, we furthermore assumed that the particles were compact, 
and the turbulence parameter (a) was set to zero. Later on we investigated the effects 
of different porosity models, more realistic collision models (with sticking, bouncing and 
fragmentation) and turbulence of different strengths. Below we summarize our results. 



• Porosity helps to produce heavier particles, and it can somewhat delay sedimen- 
tation (e.g. increase the sedimentation timescale), although this delay is limited. 
Without the stirring effect of turbulence, particles inevitably settle down to the 
midplane. 



Upon using the porosity model of lOrmel et al.l (120071 ) . the enlarge ment parameter 



of the particles is generally between 10 and 100. Upon using the 
( 2009 ) . the enlargement parameter is between 10'^ and 10^. 



Okuzumi et al 



• Bouncing prevents particles to reach masses greater than ~ 1 g (the exact value 
depends on the disk, porosity, and collision models). 

• As bouncing results in a narrow size distribution and halts particle growth, turbu- 
lent mixing equalizes the particle properties. The mass and porosity is constant as 
a function of height above the midplane, if the turbulence parameter is constant. 
Unless a dead-zone is present at the midplane, the particle properties observed in 
the disk atmosphere directly reflect those at the midplane. 

• A higher value of turbulence decreases the particle masses but it increases the dust 
scale height. Using the simplified Braunschweig model with a = results in 
a dust scale height of 0.1 Hg and final particle mass of 10"^ g, however the dust 
scale height is 0.5Hg, and the final particle mass is 10~^ g when using a = 10~^. 
Therefore, a sufficiently high turbulence value can keep the disk atmosphere dusty 
but the absence of significant dust growth is not favorable for planet formation. 

• When using the most detailed collision model up to date (the Braunschweig colli- 
sion model), we obtain particle masses of 10^2 g (^iti^ 

an average radius of 1 mm, 

and an average Stokes number of 2 x 10~^) and a dust scale height of 0.25 Hg in 
the considered disk model. 
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Conclusions and Outlook 



In this thesis, numerical Monte Carlo simulations of dust growth were performed in 
local box (OD) and vertical column (ID) models using a laboratory experiment-based 
collision model to better constrain the initial stages of planet formation. Previous dust 
coagulation simulations usually used a simplified collision model with sticking and frag- 
mentation. However, the laboratory experiments performed during the last 15-20 years 
revealed a variety of different collision types (penetration, erosion, bouncing, etc.). The 
main idea behind this thesis was to construct a collision model using all the available 
laboratory data for silicate particles and implement this collision model into a dust 
evolution code. 

We identified nine different collision types based on 19 experiments and constructed 
eight collision regimes based on the mass ratio and the porosity of the particles (similar 
sized versus different sized particles, fluffy versus compact particles). The experiments 
showed that sticking collisions mostly happen at low collision energies. At intermediate 
energies, bouncing is dominant; and if the relative velocity of the two colliding aggre- 
gates is higher than 1 m/s, fragmentation occurs. Naturally the whole parameter space 
cannot be covered by experiments, therefore we had to extrapolate at areas where no 
experiments were performed. Upcoming new experiments have to confirm whether these 
extrapolations are correct or the collision model needs to be modified. 

Calculations using previous collision models (based on sticking and fragmentation only) 
showed that an equilibrium between fragmentation and sticking is reached. In our 
collision model however this docs not happen since the particles enter the bouncing 

/)jSU9p 550LU 'LUJOU 

regime and at that point their growth is halted. This result has strong consequences. 'a o o o 
1.) As small particles arc not produced in a continuous manner by fragmentation, 5 3 2 5 S 
it is not clear what mechanism can keep the disk atmospheres dusty throughout the 
lifetime of the disk as it is observed. 2.) As the growth of the particles is halted 
already at intermediate collision energies, the particles produced in our simulations are 
generally smaller than in previous simulations of dust growth, which is not favorable for 
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planetesimal formation. We performed simulations of three different colhsion models, 
using three different turbulence strengths and we found that the Stokes number of these 
particles is rather insensitive to the disk parameters. It is always around St = 10~^. 
It has to be investigat ed whether these small particles are la rge enough for particle 



clumping mechanisms (jCuzzi et al. 



2008 



Johansen et al 



20071 ) to form planetesimals. 



Previous models of the vertical structure of disks suggested that in order to keep the 
dust atmospheres dusty for the observed lifetime of the disk (some 10^ yrs), some kind 
of grain retention mechanism is needed. However, as discussed earlier, bouncing halts 
particle growth therefore small particles are not produced by fragmentation. For this 
reason we revisited the problem of vertical settling using the new collision model and we 
were searching for ways to delay particle sedimentation. The effects of different porosity 
models were investigated and we found that porosity alone cannot sufficiently delay sed- 
imentation. A simplified Braunschweig collision model was constructed using sticking, 
bouncing and fragmentation only. Such a collision model slows down sedimentation as 
particles are kept small (10~^ - 10~^ g) and even an intermediate level of turbulence 
can stir up and keep these particles levitated. Upon using the complete Braunschweig 
model, we found that particles can reach higher masses (10-2 g) and results, the 
scale height of the dust is smaller. 

There are various opportunities for future work, most of which are highlighted in the 
discussion sections of the various chapters. Here I mention a few of them which I believe 
are of major importance. 



Planetesimal formation by coagulation: iBrauer et al.l (l2008bl ) showed that a pres- 
sure bump stops the radial drift of the particles through the disk and in these areas the 
dust particles are accumulated, therefore they can form planetesimals by coagulation. 
This problem should be revisited using the collision model described in this thesis as 
bouncing might prevent the formation of these planetesimals. 

The radial drift and stopping time of aggregates: The radial drift of aggregates 
through the gas disk is one of the biggest problems of planet formation as the solid ma- 
terial can be rapidly lost to the star. It is therefore astonishing that the most commonly 
used formulae to calculate the stopping time of aggregates (the coupling strength of the 
dust to the gas, which de termines the ra dial drift) were de rived in t h e '70 and at earlier 



times (Epstein regime - lEpsteinl (|1924| ): Stokes regime - IWhippld (|l972l )). The basic 



assumption that enters these models is that the aggre gates are spheres. Sinc e then, not 



so much work was spent on investigating the problem (IMeakin &: Donn 



1988 



_Nakamura 

& Hidaka, Il998l ). It would be worthwhile to study by numerical simulations how the 
stopping time is altered for e.g. fractal or fluffy aggregates. Using such non-uniform 
or non-spherical structures, a significant amount of the collision energy between a gas 
atom/molecule and the aggregate might be converted into rotation energy instead of ki- 
netic energy. Such calculations might show that the radial drift is reduced for aggregates 
and new, more realistic formulae for the stopping time could be derived. 
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New hit & stick poros ity model: We inv estigated two porosity models in Chapter 
[5j The porosity model of 



OrmeletaL 



(j2007l ) probably underestimates the fluffiness of 
the aggregates by not tak ing into account collisi ons between different sized aggregate 
structures. The model of lOkuzumi et alJ (|2009l ) treats such collisions, however their 
model is based on geometrical sticking (no restructuring). Therefore their model might 
overestimate the aggregate porosity. We therefore propose the development of a new 
porosity model that takes into account collisions between different sized aggregates, and 
in the same time does not assume geometrical sticking. As the porosity determines the 
coupling strength between the dust and the gas, better understanding of the porosity is 
crucial for coagulation calculations. 

Properties of icy and organic monomers: In this thesis we considered the coagu- 
lation of silicates with 1.5 micron diameter, however other materials and sizes should be 
considered, too. Performing laboratory experiments with such materials is a challenge, 
but their importance for planet formation is potentially great. Such monomers might 
have much better sticking properties than silicates, therefore these materials can turn 
out to be favorable for planet formation. Although some experiments were performed 
on macroscopic bodies (see Sec. II. 3|) . we still do not know the properties of microscopic 
bodies. 
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